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STABILITY OF LARGE SCALE POWER SYSTEMS 

Transient stability of a power system is an 
attribute of the system which denotes the conditions 
in which the various synchronous machines of the system 
remain in S5nichronism with each other when sudden 
disturbances occur on the system. Stability studies 
are necessary for planning new facilities for future 
load growth and also for a reliable operation of the 
system. For such studies it is now acknowledged that 
Lyapunov method yields satisfactory results as compared 
with the conventional method using repeated integration 
of the system eq.uations for different assiuaed clearing 
times. Ever since the original works of El— Abiad and 
Nagappan, Gless and Aylett in this area there has been 
considerable interest in the use of Lyapunov functions 
for assessing power system stability. However, after 
1972, the interest somewhat declined due to some of the 
practical difficulties encoimtered in the successful 
application of this method to realistic power systems. 
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Some of these are: 

(i) ■ Determination of the stability boundary for 

the post-fault system. This involved the 
computation of 2 -1 unstablo equilibrium 

points. This was a formidable task even for 
off line studies, 

(ii) The larger size of the system necessitated 

the use of simpler models for the synchronous 
machines, 

and (iii) Increased interconnection of the system required 
the construction of the dynamic equivalents and 
their incorporation in Lyapunov methods. 

In this thesis solution to some of these problems 
is presented. As regards the determination of the stability- 
boundary the recent work of Prabhakara and Bl-Abiad marks 
a big step forward. However, it has certain drawbacks in 
terms of choice of reference menhine and an adequate 
theoretical basis. In this thesis the problem of 
determining the stability regions as formulated in the 
context of a multinonlinear lure '-Popov system where the 
nonlinearities satisfy the sector condition only in a 
region arcund the origin of the state space. Utilizing 
this sector violation information of the nonlinearities, 
a generalization of an earlier resiilt due to Walker and 
Me ' Clamroch and Weissenberger is obtained. This result 



sii 


is further suita.bly modified to suit its application to 
the power system model. The method of Prabhakara and 
Bl-Abiad is shown to be a special case of the general 
result obtained. 

The problem of handling large scale power systems 
is analyzed by suggesting a method of decomposition and 
the use of a new theoretical concept of vector Lyapunov 
functions. The large scale power system is decomposed 
into low order models called the 'subsystems’ and the 
stability properties of these subsystems are analyzed 
separately by constructing scalar Lyapunov functions 
for each of these subsystems. On the higher hierarchia.1 
level these scalar functions a.re used to define a vector 
Lyfipunov fmction for the composite system. Using this 
vector function .together with the na.ture of inter- 
actions a.mong the subsystems, the stability of the overall 
system is inferred. A 3 machine example illustrates the 
application of the proposed method. 

A third contribution of the thesis is in the area 
of using energy or energy type Lyapunov functions to 
identify coherency and also subseq.uently develop dynairiic 
equivalents for large power systems. Stability analysis 
by both conventional and Lyapunov based methods is 
carried out for the reduced and the unreduced systems, 
and the results are compared. 



The following is a chaptoni-rise suimna-ry of the 
work reported in this thesis • 

The first chapter introduces the problem of 
stability analysis of large power systems using the 
second method of Lyapunov. The state of the art is 
briefly reviewed and then the scope and objective of 
the thesis are outlined, * 

The second chapter is devoted to the development 
of state space models of power systems. Although a proper 
rationale involving the control theoretic concepts of 
controllability and observability has been given for a 
proper choice of the state variables, the reasons for 
such uncontrollability are not evident. In this chapter 
it is shown that this state uncontrollability is partially 
due to the state overdescription in the model and hence 
the corresponding state vector is not of the proper 
order. A legitimate order of the state vector is obtain- 
ed and is shown to be the same as that obtained by the 
control theoretic approach. Thus some additional light 
is thrown on this already settled controversy on the 
minimal order of the system. 

The third chapter introduces a new concept of 
vector Lyapunov functions for the stability of a large 
scale power systems. The theory of vector Lyapunov 
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function is developed to suit, the pouor systeiri problem. 

A decomposition of the power systora into lower order 
subsystems is obtained. Each of these subsystems is 
shoTO to be identica.1 in form and in order. Lyapunov 
fuactions for each of those subsystems are constructed 
using the Moore-Anderson’ s tiioorem. These functions 
have a positive definite quadratic form with p. negative 
definite derivative. A vector Lyapunov function is 
defined for the composite system. Conditions for the 
stability of the overall system aro then derived. A 
scalar Lyapunov function is also defined for the 
G ompo si iB system and is proved to possess a negative 
definite derivative. This sca,lar function is sub- 
sequently used for determining the stability domain. 
Illustrative o.xample of a three machine case is given. 

The fourth chapter deals with the problem of 
obtaining stability region for a multimachine power 
system without computing the unstable equilibrium points. 
The problem is posed first as a multinonlinear problem 
of the Lure ’-Popov type in v/hich the nonlinearities 
violate their respective Popov sectors . Using this 
sector violation information, an expression for the 
stability boundary is obtained constituting a generali- 
zation of the results of Wa-lker o.nd Mc’Clamroch obta.ined 
for a system with a single nonlinearity. The technique 
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is tiieii extended to a nultima chine power system for 
estimating the stability domain. Two algorithms are 
proposed for practical implementations. The methods 
are illustrated with reference to a 5-ma.chine exa.mple. 

The fifth chapter a-ddresses itself to the problem 
of stability analysis using , reduced order models of 
the power system. A criterion using energy functions 
for identifying coherent group and for the development 
of the dynamic equivalent is proposed. Test results 
on a 44 bus system with 15 generators are presented. 

A general review of the results of this research, 
identification of some unsolved problems and suggestions 
for future work form the theme of the concluding sixth 
chapter . 



CHAPTER I 


DTTRODUCTIOH 

1.1 THE STABILITY PROBLEM: 

The term ’Stability’ when used with reference to 
a power system is that property of the system or a part 
of the system which enables it to develop restoring forces 
between its elements equal to or greater than the distur- 
bing forces so as to restore a state of equilibria between 
its elements [l]. (imerican Standards Definitions of 
Electrical terms, ASA-C42-1941) . Although stability is a 
single phenomenon, there are two important stability concepts, 
namely, the local stability and the transient stability. 
Whereas the former, usually referred to as 'Steady state 
stability* , deals with the stability of the equilibrium 
state of the system for small deviations of the system 
variables such as the load disturbances or prime mover 
inputs, the latter is concerned x-rith the stability and 
transmittal of power when the system is subjected to a severe 
'aperiodic disturbance*. By an aperiodic disturbance, it 
is meant that it does not come with regularity but occurs 
after an interval during which the system reaches an 
equilibrium state between disturbances [1]« These distur- 
bances may be due to sudden loss of load, loss of excitation, 
switching operation or sudden short circuits. Of these the 
short circuit is the most severe and, unless it is cleared 
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quickly, results in a loss of synchronism of the system. 

It is therefore evident that there exists a maximum time 
beyond which the failLt cannot be allowed to persist on 
the system. This is called the critical clearing time. 

An assessment of this critical clearing time is extremely 
important in planning studies, relay coordination studies etc. 

The problem of transient stability is therefore 
formulated in the following manner: Given a system, does 
there exist an equilibrium state of the system after the 
disturbance is cleared ? If yes, what is the maximum time 
that the disturbance may be allowed to remain without the 
system losing synchronism ? An analysis of the problon 
essentially involves the following two main steps; 

Step 1: The study of the evolution of the system from 
the occurrence of the disturbance to the time 
of clearing. 

Step 2; The study of the evolution of the system beyond 
the time of clearing. 

The intervals of time corresponding to these two steps are 
called the 'faulted state’ and the 'post-fault state' 
respectively. Bach of these steps involves the integration 
of a large number of nonlinear differential equations. 

The topic of transient stability investigation dates 
back to the early 1920 ’s when the point-by-point method [2] 
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was -used for the analysis of simple 2 or 3 machine systems. 
Subsequently 5 tho network analyzer was used to simulate 
larger sized systems for carrying out these investigations. 
With the increase in size and complexity of the power system, 
the digital computer has attained undisputed supremacy over 
the network analyzer for such and other studies. 

The classical method of stability investigation 
consists in assuming a clearing time 't ' arbitrarily and 
solving the differential equations of the system in its 
faulted state upto t and the post fault state after t to 
obtain the swing curves . The initial conditions are obtained 
from a prefault load flow data. These swing curves which are 
a plot of the rotor angles vers-us time are examined for the 
stability of tho post-fault system. This is usually done by 
observing whether the rotor angle differences with respect 
to any one machine taken as a reference tend to remain 
constant during the first swing of the rotor angles. In 
such a case we conclude that the system is stable. If they 
tend to diverge, then it is considered to be unstable. 
Generally it is assumed that if the system is stable during 
the first swing, it is also stable thereafter This process 
is repeated for another clearing time. After a certain 
number of such repea,t procedures, the clearing time at which 
the transition from the sta,ble to the unstable state of 
the system takes place is established. This gives the 
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critical clearing tine. Standard nmerice.1 integration 
and associated prograns are available to c'^jrry out such 
investigations. Hevorthelc-ss , this procedure of direct 
determination of tho transients is not an efficient way 
from the computational point of view when ono looks for a 
quick information about the stability with respect to 
various possible perturbations. It is because of this 
reason, specific care has been devoted in the last few 
years to examine alternatives such as direct analytical 
methods based on Lyapunov theory. 

Stability investigations by direct methods in con- 
trast to the classical methods are carried out by making 
partial or no use of the system differential equations. 

The equal area criterion [2] is a graphical procedure for 
single and two machine systems and does not make use of the 
differential equa.tions as such. Analytical methods such as 
tho phase plane technique [3] were employed for single and 
two machine systems. The energy integral criterion [4] 
and the direct method of Lyapimov [5-5 5] for multimachine 
systems make use of the system differential equations 
partially. Of these, the Lyapunov’s direct method has been 
a topic of active research in recent years. The technique 
essentially consists in establishing a region of stability 
around the post-fault equilibri-um state via the construction 
of a Lyapunov function. The system is considered stable if 
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the solution tro.j ectories of the faulted and the post-fault 
states lie entirely in this region [S]. 

1.2 lYiiJOTOY’S DIESGT ivlETHOh APPLIED TO POVJBE SYSTEMS: 

The application of Lyapunov’s direct method to the 
transient stability studies of nultimachinc power systems 
consists in the replwicenent of the step 2 detailed above by 
a stability criterion in the form -of a Lyapunov function 
around the post-fault equilibrium point. This construction 
of a Lyapunov function, in fact, is a key step in the appli- 
cation of the Lyap'unov theory to the power system stability 
problem. Such a criterion was first suggested by Aylett [4] 
in 1958 based on energetical considerations but the theory 
was not well developed then for practical implementation. 

The pioneering work of Bl-Abiad and Nagappan [6] and G-less 
Gless [5] triggered a. series of investigations in this area. 
An excellent snrvey of the literature in connection with 
the construction of Lyapunov function is contained in 
references [22 ]and [11 ]. The latter also contains a critical 
review of these functions. Two different linos of approach 
have been adopted by researchers for the construction of 
the Lyapunov functions. They are (i) the Ealman's 
construction procedure used by Pai et al [ 9 j and Pai [7] 
and (ii) the Anderson's technique applied by Willems [14,15 ] 

and Pai and Murthy [lO]. In reference[lO ],a Lyapunov 

T 

function was obtained with a quadratic part X PX which 
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lias a positive definite syanetric P. Such a function 
has "been found useful in the application of vector 
Lyapunov function theory [52] which constitutes a part 
of this thesis. Recent works of Mans our [l9] and 
Willens [2l] generalize the results so fa.r obtained in the 
literature via generalized Lyapunov functions. 

The procedure of applying Lyapunov's second method 
to the power system stability problem consists in obtaining 
a suitable mathematical model for the system of the type 

X = I (X) (1.1) 

for both the faulted and the post-fault states, let these 
be represented by 

i = (X) (1.2) 

and X = ^2 (2l) 

for the respective states. A Lyapunov function V(2) is 
them constructed for the system (1.3). Then a stability 
region around the origin of the post-fault ociuilibrium state 
X = 0 is defined by an ine equality 

V (X) < e (1.4) 

where e is a constant to be computed. This constant 
is usually obtained by evaluating the function V (X) at 
the unstable equilibrium point closest to the post-fault 
stable equilibrium point. Now the calculation of the 
critical clearing tine is simple and is done by integrating 
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the differential equations (1.2) and ev.al’aating '/ (X) at 
every tine step of the integration. The tine t at which 

7 (X) = £ (1.5) 

gives an estinate of the critical clearing tine. Ihus an 
explicit integration of the differential equations (1.2) 
is perforaed only once imlike the classico.1 techniques. 

This results in a, significant reduction in the cenputation 
tine particularly in the analysis of large sca.le power 
systems for various contingencies. In this case one can 
consider the Lyapunov method as complementing the classical 
methods. A preliminary screening by this method can reduce 
the number of studies to be carried out in dota.il. Another 
advanta,ge of the method lies in investigating the effects 
of par.aaeter variations such as damping [ 7j 14) » power 
inputs [ 23 ] etc. , thoroughly and quickly. Stability 
indices can also bo defined by this methcc: without making 
actual stability calcu.laticns . It is now groxlually being 
recognized that the method is able to give satisfactory 
results in agreement with those obtained by dix'oet simulation. 

1.3 THE LARGE SCALE POVER SYSTEM STiiBILITY PROBLEM; 

The technique of applying Lyapunov methods has 
been described in the previoois section. Maile consider- 
a,ble efforts have been directed towards the improvement 
of the Lyapunov functions [10,19, 21] on the one hand, the 
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application of tho tliecry to realistic systcns on tiie 
ether hand has been very slew inspite of the eidvantagos 
mentioned earlier. This is due to the fact that there 
are certain drawbacks, apart from the general conservative- 
ness of the method, thiat still need further a,ttention. 

Firstly, the application of the theory to largo power 
systems poses heavy computational problems in the deter- 
mination of stability regions [35]. It has been established 
that for an n-machine system, the number of unstable 
equilibrium points is 2^“^-l, and consequently very large 
for a large 'n'. Although efficient algorithms are 
available [12,13,26,28], considerable computation time is 
required. This problem has been solved to a great extent 
by Prabhakara and El-Abiad [31] recently but needs further 
refinements tc yield better stability regions. Another 
drawback of the method is in regard to the modelling of 
the system. The selection of the models for the system in 
Lyapunov methods is strongly limited by tho practical 
possibility of constructing a suitable Lyapunov function 135 ]• 
Unfortunately, till now it has not been possible to 
include refinements like governor dynamics, saturation 
and saliency effects etc,, into the Lyapunov functions 
for multimachine systems. Therefore the general problem 
of multimachine systems has been considered [13] only on 
the basis of certain simplifying assumptions, Tho 
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scluticn tc scriG of these fcregcing probleris calls for 
a new approach tc the prchlen ...f stability cf large scale 
power systems. In this thesis an a.ttenpt has been made to 
cover some new grounds in this direction. They are 

(i) Application of the concept cf voctor 
lyaptinov function for the large scale 
power system stability problem. 

(ii) Development of new algorithms tc doter- 
mine the stability boundaries of multi- 
machine systems 

and(iii) development of d3naamic equivalents for 

power systems xasing the Lyapunov function 
concept. 

¥e now briefly review the pertinent literature in these 
areas before outlining the objective and the scope of the 
present thesis. 

1.4 THE STATE OP TUB ART: 

1.4.1 The vector Lyapunov function approach: 

Most complex systems are made up of a number cf 
low order interacting systems. These low order systems 
in hierarchial theory [ 36J are called 'Subsystems’. These 
subsystems are arrived at on the basis of docompositicn of 
the large system. Decomposition and aggregation have 
long been used in economic theory and found an engineering 
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application [37] tlircugh the tearing techniq_ucs of 
Ercn [33] • Today's power systens being largo and ccnplex 
becoEe natural candidates for such a.pplicp.tion. The 
analysis of load flow and short circuit prcblens in large 
scale systens -using the diakoptic nethod is well kno-wn. 

One therefore is notivatod to extend these ideas to the 
stability prcblens of large scale systens using lyapuncv 
netheds. Such a procedure vraB first presented by 
Bailey [39]. The tochniq.ue consists in deconposing the 
systen into low order subsystens. Suitable lyap-uncv 
f-unctiens are constructed in the Tjsual nanner for these 
subsystens satisfying certain sign definite properties. 

On a higher hiera-rchial level, the concept of a vector 
Lyapunov function [40] is introduced in which these scalar 
f-unctions fora the elonents of the vector lyapimov function 
defined for the conposite systen. Using the properties of 
interactions anong the subsystens conditions for the 
stability of the overall systen are derived. Those ideas, 
however, did not receive s-ufficient attention "until recently 
when Siljak and others[41r5l] developed the theory f-urther 
in a nanner suitable for a,pplication to physical systens. 

In spite of those developnonts applications to realistic 
systens are yet unkno-wn in the literature in general and 
the area of the power systens in particular. Probably the 
reason for this is that a suitable deconpositicn procedure. 
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which, is a key concept in hierarchia.1 theory, cciild not 
be developed for inpleoentaticn cf the thecry. In this 
thesis an attenpt h?;.s been made tc apply voctcr lyapmov 
fmaction to the stability an.alysis of large scale power 
systens [52]. 

1*4.2 Stability regions: 

The problen of estino.ting the stability regions of 
nultinachine power systens for iyapunov nethods was identified 
with the conputation of the tuastable equilibrim points of 
the post-fa^nlt systen [5,6,12-25,26—28, 31] and eva,luating 

the V-f-unction a,t these points. The nininm of the values 
sc con put ed yields a constant e that gives an estinate of 
the stability dcnain defined by the' inequality (1.4). It 
is well known that there are (2’^“^-l) unstable equilibriun 
points for an n-nachine systen and the conputation of these 
points becones prohibitive when 'n' is l.?xge. In this 
connection special studies have been co,rried out in this 
area and also psoudo randon search techniques have been 
proposed. Notable onong these are the works of Ribbens 
Pavella [12], Tavcra and Snith [26,27] and Uyenura 
et al [28]. Sone interesting results based on physical 
energy type approxinate considerations [ 17,20] have also 
been proposed and applied to estinate the stability domains. 
An approxinate ’n-donensional cube* method has also been 
proposed and used [29]. The basic philosophy of this 
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neth-od ccnsists in- assiming an n-diDensional cube (n being 
tbe nmber of nachines) arcund the post- fault stable 
equilibriun point. The nininun value cf 7 (X) on this 
surface is naxinized with respect to the side length 
itself and the resulting value is used as an estimate of e 
in equation (1.4). L further approximate nothod [50,54] is 
based on a series expansion arcund suitable known points. 

The computation of V (X) at the corners of an (n-l) 
dimensional parallelopiped and choosing the minimum of 
these values was suggested by Murthy [24]. Recent work of 
Prabhakara and El-Abiad [51] approximates the unstable 
equilibrium points on considerations of those similar to 
a single machine connected to an infinite bus system. The 
present thesis utilizes the information associated with the 
nonlinearities occuxing in the mud-tinachine power system 
models cast in the lure'-Popcv form. These nonlinearities 
are known tc viola.te the Popov sector conditions. Utilizing 
this property explicit expressions are obtained for the 
stability regions based on certain approximations and 
thus constitutes a generalization of a similar work due 
to Pai et al [9] for single machine power systems. Results 
of reference[5l]are shown to be a special case of this 
general approach. 
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1,4.5 D 3 manic Equivalents: 

In this section we review the state of the art in 
the area of systen sinplification. The growth of power 
systems in the recent years has made stability analysis 
of the systen in a detailed manner difficult duo to the 
limited core memory of even the present day computers. A. 
possible solution to this problem is to reduce the size of 
the systen by the use of equivalents to represent portions 
of the network beyond the area of immediate interest. 

Usually in power systems, a fault in one section of the 
system may not significantly affect the behaviour of some 
other section of the S 3 rstem, In such cases it is not 
necessary to model this latter section in detail and can 
be replaced by an 'equivalent' representing it while studying 
the behaviour of the faulted section. Thus a reduced 
representation of the systen called the 'd 3 airmic equivalent' 
is available. One method of obtaining such an equivalent 
is by utilizing the property of ' coherency' among machines 
in the system. and replacing such machines in the equivalent 
by a single machine. Such equivalents were advocated in 
the early fifties vrhen the systems at that time were too 
big to be manageable on the A.C. network analyzers. In 
this method of equivalencing, groups of machines that 'swing' 
together, during the period of study in the faulted section, 
are combined. Identification of these coherent machines 
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is done ly conducting a detailed study cn the systen 
and applying a set of criteria [53-59] . later a dynanic 
equivalent is developed ly using certain ctlier criteria 
[53-58, 59 ]• Recognizing that the likelihood of machines 
swinging together [ 2 ] is increased hy (i) a decreased 
impedance between the machines (ii) the initial angular 
position being in close proximity (iii) the inertia constants 
being almost equa.l and (iv) remoteness of the fault or the 
source of disturbance, Lee and Schweppe [ 59 ] developed the 
'features' of pa,ttern recognition and used then for 
coherency identifica,tion. 

Other approaches in the development of dynamic 
equivalents employ modal techniques and system order 
reduction [60,61,62 - 64] • In this thesis Lyapunov functions 
are used for coherency recognition and also to obtain 
dynamic equivalents. 

1.5 SCOPE ML OUTLINE OP THE THESIS; 

The main objectives and chapterwise summary of the 
work reported in this thesis will now be outlined. 

The main objectives of this thesis are' 

(i) to present an alternative technique for the 

stability analysis of large scale power systems 
using vector Lyapunov functions, 

(ii) to derive explicit expressions for the stability 
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donains of nultiaaotiiiiG power systcns based 
on the sector violaticn criteria of nulti- 
nonlinear systens, 

and(iii) tc develop an effective criterion for 

coherency reccgnition and ccnstruction of 
a d3na.9inic equivalent based on lyapunov’s 
nethod. 

The following is a chapterwise suinnary of the work being 
reported in this thesis. 

Chapter II begins with the definitions of ovcrdescribed 
models [65,66] of systems and relates these concepts to state 
uncontrcllability. The state space model of power systens 
in the 2n dimensiono,! state sijace is. developed for a 
2-nachine power system and is shown tc be overd escribed. 
Explicit dependency relationships amcng the state variables 
are derived both for uniform and nonuniform damping cases, 
from which the legitimate or the minimal order cf tbe state 
vector is derived. The analysis is then genoreilized to 
an n-nachine system. State models in the 2n. and (n+m) 
dimensional state spaces for multimachine systens for the 
uniform and the nonuinifcrm damping cases respectively are 
then introduced. These models are useful in decomposing 
the multimachine system into lower order subsystems for 
the application of the vector Lyapunov approach in 
Chapter III. These are shorn, to be overdescribed and the 
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legitimate order of the state vector is derived. It is 
ccncltided that the legitimate order of the state vector 
for an n-machine system is (2n-2) in the case of -uniform 
damping and (2n-l) for the nonuniform damping case. Thns 
a physical insight into the reasons of -uncontrollability 
of certain power system models is provided in this chapter. 

Chapter III introduces the concept of vector lyapunov 
function as applied to a large scale system stability problem. 
Firstly the concepts of the subsystems and their interactions 
are introduced by a state variable decomposition of the 
original dynamical system. The required conditions on the 
scalar f-unctions for asymptotic stability of the individual 
•free’ subsystems and. also those on the interactions are 
next reviewed. A vector lyapunov function for the composite 
system is then defined and conditions for the asymptotic 
stability of the overall system are given by means of a 
theorem due to Grujic and Siljak [48], This vector lyapunov 
f-unction is then used to construct a scalar Lyapunov function 
for the composite system. This latter function is usef-ul 
in determining a region of stability for the overall system. 
The application of this theory to the power system problem 
is next considered. The decomposition procedure of the 
system is outlined for a ^-machine case first and is later 
generalized to an n-^nachine system. It is shown that 
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each of these subsysteas is described in its- nininal state 
space, lure^ type lyapimcv fimcticns are constructed for 
each of then. Conditions, for stability of the post-faiilt 
composite power system are later derived via. a, vector 
Lyapunov function, A scalar Lyapunov function for the 
composite system is then constructed using the vector 
Lyapunov function to estimate the stability region. The 
entire procedure is finally illustrated with reference to 
a 3-Eiachine example . 

Chapter IV is devoted to the estimation of stability 
domains for multinonlinear systems cast in the Lure' -Popov 
form in which the ncnlinearities violate the Popov’s sector 
condition. Using this infoimation an expression is derived 
for the stability domains constituting a generalization of 
the result for a system with a single nonlinearity [ 67,68], 

An application of this result to the power system stability 
problem is next discussed. Per this purpose, the problem 
of obtaining the stability regions for multimachine system 
is formulated based on the sector violation properties of 
the power system nonlinearities. The method of application 
of the above mentioned result is then described, iin explicit 
expression for the stability region is also obtained. It is 
also shown in this chapter that a method recently described 
by Prabhakara and Sl-Abiad [ 31 ] for estimating the unstable 
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equilibri-un point can be brought into the franework of 
the present analysis, thus providing a theoretical basis 
tc their nethcd. An oltorna.te procedure tha,t provides 
better, stability regions is also proposed in this chapter. 
Finally these nethods are inplenented on a 5-nachine system 
end the critical clearing times obtained by these methods 
are compared with that determined by direct simulation. 

The application of Lyapunov functions for coherency 
identification in a multimachine power system is demonstrated 
in chapter V. Using the kinetic and the potential energy 
components of an energy type Lyapunov function a new 
definition bf coherency is proposed. This definition is 
them used for both coherency identification and development 
of a dynamic equivalent. An algorithm for coherency 
recognition on the basis of the new definition is presented. 
Later, explicit expressions for the self admittance of the 
equivalent generator replacing a coherent group and the 
transfer admittances between this equivalent machine and 
the rest of the machines are derived. The coherency 
recognition algorithm is implemented on a 44 bus 15 generator 
system. After obtaining the coherent groups," the dynamic 
equivalent is developed. Swing curves of the retained 
machines are obtained by conducting a stability study 
using this equivalent and are compared with those obtained 



from a detailed study. Use of these dynamic 0q.uivalents 
in the study of large scale systems via Lyapunov methods . 
is discussed. 

The concluding sixth chapter highlights the major 
contributions of this thesis and certain problems encouhtered 
during the course of this research are presented. Suggestions 
for future work are then given. 



CHAPTER II 


STATE OVERDESGRIPTION. ARE IMCOHTROLLABIIITT 
IIT 20mR SYSTEM MODELS 

2.1 IHTRODUCTIOH; 

This chapter deals with the development of state 
space models for multimachine power systems for use in 
transient stability analysis. A number of state models 
have been cited in the literature [ 7,10,12,14,15,16] with 
the result that there was a great deal of controversy 
regarding the legitimate order of the system to be used in 
the systematic construction of lyapimov functions. This 
controversy has, however, been settled either through the 
control theoretic concept of minimal realization [69 ] or 
by physical argnmonts as suggested by Ribbens Pavella [ 11,1^ 
It is by now well established that for constant mechanical 
inputs and constant flux linkages, a uniformly damped and 
an undamped n-machine system should be treated in a (2n-*2) 
dimensional state space and the non-uniform damping case 
of the system in the (2n-l) dimensional sta^te space. The 
question of studying all these cases in a unified frame- 
work for the sake of comparison and convenienco was analysed 
in depth by Willems [ 21] hy defining stability not as a 
state stability concept but rather an output stability 
property or what is called partial stability. 

Althou^ certain criteria ezi^t for testing the 
controllability of the system, little insight into the 
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exact mechanism which causes this uncontrollahility in 
certain of the state models is provided hy these techniques. 
There, in fact, appears to be a lacuna in justifying the 
choice of the state variables while obtaining the state 
description of the power system differential equations. The 
presence of uncontrollability is normally associated with 
the inefficiency in the manner in which the effect of the 
control propogates through the system [-TO]'. ¥hen the swing 
equations are put in the form of a set of first order 
differential equ8,tions , the power angle nonlinearities 
constitute the 'Control variables', in the language of 
modern control theory. Controllability as a physical concept 
is therefore not meaningful although it has been applied in 
a mathematical sense, rather mechanically [69]. In this 
chapter we take a view point following Johnson [65] in 
connection with ovcrdescribed system models. It will be 
shown that some of the power system models are, in fact, 
overdescribed resulting in a functional dependency of some 
of the state variables on others. It is this fa.ct that 
results in the mathematical 'Uncontrollability* of the system. 
The process of obtaining a legitimate order of the state 
vector is also described and this order is found to be 
identical to that obtained in the literature [ 69 ]. Such 
a realization is indeed a prerequisite for applying the 
results of modem stability theory. 
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2.2 IMCONTROLLABILITY AND OVERDESCRIPTION: 

The notion of ccntrollahility of a linear dynamical 
system has been "used eztensively in the control literature 
and is a. property associated with the effect of control 
forces on the state of the system. The following definition 
of controllability is due to Kalman [ 71] . 

DEEINITION 1: A system is said to b e controllable if it 

is possible to find a control vector iXt) which in a ' 
specified time t^ will transfer the system between two 
arbitrarily specified finite states and 

Clearly, a, system is said to be uncontrollable if it 
is not controllable. It is necessary to ascertain this 
property of the linear dynamical system before one can apply 
optimal control theory. However if the system is uncontroll- 
able, reasons for such unccnt reliability have not been fully 
explained. Recent work of Johnson [ 65] explains a particular 
mode of uncontrollo.bility arising due to 'System Over- 
description', ¥e thus have the following definition. 

DEFINITION 2: A system described by 

i = 1 (2.1) 

where X is an N-dimensional state vector and U is an 
r-dimensional control vector, is said to be overdescribed 
if its state vector of order N has components which are 
functionally dependent and N xp where 'p' is the legitimate 
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order of the state "vector . 

DEEIFITION 3: A state vector X is ccnsiderod to be of 

legitimate order if it conta,ins the minimal nmber of state 
variables which are functionally independent to describe 
the physical system. 

Definition (2) implies that there is a redundancy in 
the state variable description which may not be explicitly 
exhibited in the system (2,l), In such a case there always 
exists a set of functionally independent algebraic fmctions 
{ }» i=l,2,..(I-p) such that 

, ,S]j) = constant , j=l, 2, . , (N-p ) . (2.2) 

Althou^ these equations are not known a priori, the effect 
of these is manifested in (2,1), Hence the right h,and sides 
of (2.1) shculd satisfy the following set of relations: 

(Pll(X) P3_(X,U) + •• 

+ - 0 ^ 

5 (2,3) 

5 

• • • • , • jj 

where the P^'s are the components of P(X,IJ), This set of 
equations is derived by differentiating (2.2) which gives 
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N -id)* 

^ 3 X ~ ^ > 2 , . . . (N— p ) • 

i=l i 


Define 



Tlien with (2.1) and (2.5), 



931(2) • 

eq_uations (2.4) yield 
= 0 , j=l,2 ,,..(h-p) 


(2.4) 

(2.5) 


( 2 . 6 ) 


Hence (2.3) follows from (2.6), 

Again the qj-.'s o-re not known a priori. These 
qj . . (X) fom a set of Pfaffian forms of the type 


coefficients 


q)^.l(X)dXi + + • • * + 

j=l,2 ,...(H-p). (2.7) 


If these are integrahle, then, the required functionally 
dependent relations among the state variables are obtained. 
¥e are thus led to the following theorem. 

THEOEEM 2.1: [65] 

The ma.thena,tica.l model (2.1) is uncontrollable 
if there exists a positive number p<E and a set of 
coefficients 9 . . (X) , j=l, 2, . . . (N-p) , and i=l,2,.,.K, not 
all zero, such that, the following conditions are 
satisfiedi 

(i) the right hand sides of (2.1) satisfy (2.3) 
for all X and for all U 

and (ii) the coefficients 9 .. (X) form a set of Pfaffian 
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foms that are integrahle. 

In that case the uricontrollahility is due to the state 
overdescription in the sense of definition (2), 

A NOTE ON PEAEPLIN POEMS: 

A single expression of the type 
9 l 2 _(Z)dli + (p^ 2 (Z)dX 2 + . . + (p^(X)d2^ = 0 (2.8a) 

is called a Pfaffian fom. -The Pfaffian fom is said 
to be integrable if and only if there exists a scalar 
function 

5 (X) = constant (2.8b) 

such that the total derivative d? is proportional to 
(2.8a). In other words 

= 'n:(X)<p 3 _j_(X) (2.8c) 

where •n:(X) is a ccncion factor of loroportionality. Por 
cases where constants, (2.8a) is always 

integrable [65]. 

2.3 MilTHEMilTIGAI MOBEl OP THE POWER SYSTEM': 

The following simplifying assumptions are made 
in the transient sta,bility of power systems: 

(i) The voltage behind the transient reactance is 
assTMned constant. In other words it is assumed 
that the flux linkages of the various machines 
are constant. This assumption idi generally valid 
for the first swing stability analysis. 
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(ii) Dciiiip ing is o^sunsd "to bs dijcecily px'opoi'iional 
to the slixD velocity and is thus mainly due to 
the nechanical friction and aS5mchronous torques . 

(iii) The nechanical power inputs to the synchronous 

machines ore assumed constant during the transiont 
period. This is a satisfactory assumption beca'use 
the time constants of the governors a.re much 
larger than the duration of the transient swings, 

(iv) Transfer conductances are neglected. 

(v) Saturation is neglected. 

(vi) Only round rotor machines are considered thias 
neglecting the effect of saliency. Machine 
resistances are neglected. 


Under the above assumptions the differential 
equations that described the dynamics of an n-machine power 
system are given by 

d6^ 

dt^ **" ~ ~ ^ei * i=lj2,,.n. (2.9) 


The electrical power output of machine i is given by 

^ei ^ ®i ^ii “ sin(6^ - 6j) . (2.10) 


Substitution of (2a0) into (2.9) yields 
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M. 


dt 



a6, 

=1 -dV = 

n 

-.IS 

,2 * 

*1 = 

^ni - 4 hi 

• 


'j id 


( 2 . 12 ) 


The eciuilihrim states of the system (2.10) are chtained 
hy setting both 6^ and 6^ equal to zero resulting in 

= 1 % *^0 = i 3 

Althou^ there are n equations in (2.13) 

lar difference variables these equations are not over- 

determined. In fact, (2.13) suggests that 


n 


I Si = ° 

i=l 1=1 


(2.14) 


is a necessary condition for the existence of an equili- 
brium state [14]. liiis is logical sjjaee equations (2.14) 
specify the real power balance in the system, issuming 
(2.14) to hold, if all the terms in (2.13) refer to the 
post-fault system, then the post-fault rotor angle 
6?, i=l,2,...n satisfy (2.13) identically, i.e. 

p - I B. B. B... s3^(6| - = G » 

W i=l, 2,...n. 

Then equations (2.11) can be written in the fom 


(2.15) 



28 


a^6, aa,. “ 

% ^ B. . sin(6“ - a°) 

i^i- 


- ®ij - 5j) 


I \ Bj [sin(6j^ - 6^) - Biii(6? - 6p] (2.16) 


Equations (2.16) can be cast into a compact vector matrix 
differential equation of tlie form 


X = AX, - B f ( a) 0 

O' = G^X I 


(2.17) 


where the state vector X . is of the proper order depending 
on the choice of the state variables. So too, are the 
matrices A,B and G. The vector a is an m (= n(n-l)/2) 
dimensional vector and f(o) is a vector valued function 
with m components whose i'*'^ component depends on the i"^^ 
component of the outx)ut vector a. Much of the earlier 
controversy centered aro-und the dimensionality of the 
state vector X and the choice of the state variables with 
the matrices A,B and G dependent on the way the state 
variables are selected. We shall new consider the diff- 
erent forms of the state representations discussed in the 
literat'ure and verify whether or not the choice of the 
corresponding state vector results in an overdescription 
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of the dynamical system. The case of a 2~machine 
system will be considered first and later generalized 
to an n-machine system. Explicit relationship existing 
between the states of the system due to overdescription 
will be obtained which in fact results in the so called 
•Uncontrollability' of the mathematical model. It may 
be emphasized once again that overdescription implies 
a redundancy in the description of the state variables 
leading to a functional dependency among these variables. 


2.4 OVERDESGRIPTIOU, EBDUNDMCY AUD UUGOUTROLli'iBILITY: 


2.4,1 Two machine case: 

The following choice of the state variable was 
initially found to be an obvious one [14] to put the 
system description in the form (2.17) 



i' 


(2.18) 


The state equations of the 2-machine system from (2.16) 
and (2.17) with n=2 are given by 


X. 


X, 


X, 


X, 


-X 


1 

0 -X, 

I 

1 0 


0 0 
0 0 
0 0 



1 

: ' 

l/Ml 


C\J 

1 

1 

H 

ro 


X, ' 
5 


0 


I 

. 1 

A 

0 

• 


f(a) 


(2.19) 


0 


10 0 
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a = [ 0 0 1 -1] 


^2 

X, 

where f(a) = B2 ^12 + (<^1 “ 6^ )-sin(6j-62)] 

The two machine system is thus cast in the 2n (=4) dimen- 
sional state spa,ce. 

Eq^iiations (2,19) can he put in the compact form 

i = (2.20a) 


where B is a 2n -dimensional vector function, the elements 
of which are given hy the right hand sides of (2,19)* 

Thus 

I'ld) = -A^x^ - f(0) 

(2.20h) 

B^d) = X^ 
and I'4(^) = ^2 • 

If the system is overd escribed, this set of eq_uations 
should satisfy a set of simultaneous linear algebraic 
relations of the type (2.3) and given by 

9ji(^)J'l(l) + ••• -*■ = 0, 

3 = 1,2,«,( 4— p ) ( 2 , 21 ) 

where *p’'is the imknown but legitimate order of the 
state vector X. 
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Substituting (2.20) into (2.21) we obtain 

9ji(X;){-X3_X^-(l/M^)f(cr)> + (p^.2(X){-X2X2+(l/M2)f(a) } 

+ 9j3(X) X3_ + cp^.^(Z) X2 = 0, 3=1,2,... (iT-p). (2.22) 

It is found that only one set of (p-.’s, 1=1, 2,... 4 satisfying 
equation (2,22) is possible. This set is given by 

= ^2 ' 933(1) = and cp.^(X) = D2 

(2.25) 

These 9-4^’ s form an integrable Pfaffian form given by 

M3_dX^ + M2dX2 + D3_dX^ + ^2^^4 = ^ • (2.24) 

This implies 3=1 and hence p = 4-1=3 is the legitimate order. 
Integration of (2,24) yields 

M^2l + M222 + D3_X^ + D2X4 = E (2,25) 

where X is a constant, -This is a functional dependency 
relationship among the state variables of X, This in a 
geometric sense implies that, the set of all the physical 
states { X^,X2 ,X^,Z^ ^ lie on a 3-dimensional manifold M in 
a 4-dimensional state space as shown in Pigure 2.1, States 
of the system not on M do not physically exist and hence 
cannot be reached by the physical system irrespective of 
the way in which initial states are chosen. This, in the 
mathematical sense results in uncontrollability of the 
physical system. Application of the standard criteria 
for controllability to (2.19) will also reveal it to be 




33 


uncontrollable. It is therefore necessary to redefine 
the state variables. Two different choices are possible 
(i) relative rotor velocities and (ii) relative rotor 
angles. At this juncture two different cases of the 
system arise. They are 

(i) Fon-unif orm damping where ^ X 2 
and (ii) uniform damping where X^ = X 2 . 

We shall consider each of these cases separately. 

Case (i)s Non-uniform dampings 

Equations (2,9) show that the state model cannot 
be constructed with relative rotor velocities as the 
state variables thus leaving relative rotor angles to be 
considered for state representation. Thus the legitimate 
sta,te vector has components given by 



X2 - w 2 

and • 


(2.26) 


With this set of variables, the state space model is of the 
form 


il 



0 

0 ' 


— 1 

h 


1 

H 

H 

1 


ss 

0 


0 


X, 

1 

-I/M2 

•M 

1 


1 

-1 

0 



1 

0 
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a = [0 


0 


1 ] 




(2.27) 


It is easily verified in a manner outlined in the 
previous section that no nontrivial found for 

the system (2.27) satisfying (2.3). Thus the legitimate 
order is indeed 3. Also, the usual controllability 
conditions reveal that the system is completely controllable. 
The transfer function ¥jj(s) of the linear part of the system 
is given by 


^^(s+X^) + 


(2.28) 


Subscript M refers, to the nonrruniform damping case. 


Gase (ii)2 Uniform damping! 

By setting X^ =\ 2 - ^ equation (2.27) one obtains 
a system description of the form 


" V 


1 

1 

>* 

0 



C\J 

= 

0 

0 


X 2 

X^ 


0 * 

H 

1 

H 

J 


3_.: 


[0 0 1 ] 


x. 


X, 


1/M3 

-1/M. 


f (o) 


(2.29) 
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Por these equations no separate set of <p..’s can be 
obtained that satisfy (2.3). Hence the system is not 
overdescribed in the sense of Johnson [65]. But from 
(2.29), one can trivially obtain a relationship of the 
type 

+ ^2^2) (2.30) 


which can be easily integrated to give 


+ M2X2 ={ M^X3_(0) + M2Z2(0)} (2.31) 

where Z^(0) and 12(0) are the initial values of the states 
Y and X 

1 2. This exhibits a clear dependency among the 

variables. This sort of physical dependency relationship 
among the state variables was pointed out by Rosenbrock [66] 
exhibiting thus a clear redundancy in the state variables . 
Johnson’s conditions fail to show this sort of over- 
description because theorem 1 is only a necessary 
condition for overdescribed systems. It is also observed 
that this pehnomenon manifests itself in pole zero 
cancellations in the transfer function Vir^(s) of the system. 
Prom equations (2.28), by setting =X2 =X it is 

noticed that there is a pole zero cancellation. Thus 



¥y(s) = 


s( s+X) 


(2.32) 
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Subscript IT refers to the uniform damping case. 

This redundancy again calls for a redefinition of the state 
variables. The relative velocity ^ 2 } chosen 

as one of the state variables instead of and 0 ) 2 . Thus 
with relative velocity and relative rotor angle as the 
state variables, the state vector is defined by 


h = “i - "2 

Xj = - 62 ) - ( 6 ° - 62 ). (2.33a) 


Eq.uations (2.9) can be recast in the form 


^ ^ 1 ' 


-X 0 ' 


h ; 


, 1/M^+1/M2 



: 1 0 ,, 

_ j 


Xj , 
J 


0 


cr 


[0 



(2.33b) 


f(cr) is the same as in equation (2.19). It is easily 
verified that (2.33) is completely controllable and hence 
the order of the legitimate state vector in this case is 2. 


REMARK: The state models derived, in equations (2.29) and 

( 2 . 33 ) are useful in constructing Lyapunov functions for the 
subsystems in Chapter III which deals with vector Lyapunov 
functions. 


¥e now generalize these results for an n-machine 
power system described in the 2n-dimensional state space 
[7,8,14,15]. 
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2.4.2 n-macliine cases 

If we define the state variables by 


Zi =<0 . 


, i =1,2,... n. 


(2.34) 


the state model is then given by 
i = A*X - B*f(a;) 

where 


A* = 


-A 

0 



^ -1 1 

nn 

nn 

o 

9 

B* = 

nm 

I 

0 



0 

nn 

nn 



nm 








(2.35) 


and 5^ = [0^^ K^] 


in which A = diag(X^), i = l,2,,.,n, M = diag(Mj^) , 
i=l,2,...n. 0 ^^ and are the nnll and unit matrices 

of order n respectively. The matrix (m = n(n-l)/2) 

is given by [7] 

K, 


( 2 . 36 ) 


= [K. : K, ; . .. :k . ] 

nm ‘-1.2) I n-i-' 

where i = l,2,...(n-l) is an nx(n-i) matrix having the 

following structure [7] 


= 


‘^(i-1) (n-i) 


1 1 

-1 0 

0 -1 


1 

0 

0 


L 


0 0 


-1 


X2.37) 
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Double subscript notation in equations (2,35) - (2.37) 
is used to identify the rows and columns of the respective 
matrices. The function f(a) is an m-vector \fhose ith 
component is a function of the ith component of 0 only. 
Thus ^jL('^j_) is given by 

^i^^i^ = \ Ej Bj^^[sin( 0 j_ + 6^) - sin i = l, 25 ,...m 

( 2 . 38 ) 

where 9 = 6° and 6° is an n-vector given by 

6° = 6°]^. (2.39) 


It has been shown by Murtby [24] that the system (2.35) 
is uncontrollable and unobservable. Choose the 9 jj_’s 
such that the H-vector (here K = 2n) constituting the 
(p^^'s as its elements and of the form 


'P32 




T 


(2.40) 


is given by 

£ = M 2 , ... M^, D^, D 2 , ... . (2.41) 

It can now be shown that 

0 ( 2 . 42 ) 

where P is a vector of the functions on the right hand 
sides of ( 2 . 35 ). These cp^-,-’s yield an integrable Pfaffian 

J J- 


forms of the type 
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n 

I 

i=l 


f _ M.dZ. + D.dZ^^^ = 0 


On integration (2.45) results in 


(2.45) 


n 

I (M.Z. +D.Z . ) = constant (2.44) 

thereby concluding that the system is oTerdescribed in 

this state space. For the non-imif orm damping case, 

(2.44) is the only functional dependency relationship and 
hence the legitimate order 'p' of the state vector Z is 
given by 


p = (2n - 1). (2.45) 

For the uniform damping case, however, the order of 
the statevector can be reduced further by 1 in view of 
the relation 

I M-Z. = ( I M.X.(0) e"^^) (2.46) 

1=1 ^ ^ i=l ^ ^ 

which again exhibits a redundancy in the state variable 
description according to Rosenbrock [66], Hence the 
legitimate order of the state vector is (2n-2) for the 
uniform damping case. 

2.5 ALTERNATE FORMS OF n-MAGHIHB MODELS s 

We shall next consider two particular state 
models in the 2m (m = ^(n-1)) and (m+n) state spaces for 
the uniform damping and non-imiform damping cases 
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respectively. In view of the discussion in Sec. 2. 4 
regarding the legitimate system order, these models are 
overdescribed. However, their consideration is merited 
by the fact that these models lend themselves for decompo- 
sition into subsystems for applying vector lypunov functions. 
Together v^ith the description of these models we briefly 
indicate the nature of overdescription. 

2.5.2 System description in the 2m-diinensional 
state space [12]: (Uniform damping only) 


The state variables chosen to represent the 
uniformly damped case [12] in the 2m(=n(n-l) ) dim.ensional 
state space are given by 


h = “k - “o 


h-m = 


A 

5 


(i=l, 2, . . .m and 


( 62 - 6 °) 


5 k, j = 1,2, 0 ..n,k < j) 


The state space model is then given by 
Z = A 2 - B f (o) 
o = C® Z 


(2.47) 


(2.48) 


where 




— A •* 

0 


mn 

mm 

A = 

1 

^mm 

^mm 

and 0 




I B 


k'^ M“^ k ^ 

nm nm 


0 . 


mm 
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Th.8 fj^(a^)'s are the same as in equation (2.38). Over- 
description in these models will now be illustrated with, 
reference to a 3-machine example. 


3-Machine cases 

Consider a 3-raachine system with uniform damping. 
The state variables according to equation (2,.47) are 

X2 = (^1 - 0)3) .52^ = (6^-6^) - (6^-6^) 

= (a>2 ; Xg = - (62“*^3^ 

(2.49) 

since 2m is 6. The state model is given by 
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and 


^1 


0 

0 

0 

1 

0 

0 

^2 

= 

0 

0 

0 

0 

1 

0 



0 

0 

0 

0 

0 

1 


( 2 . 50 ) 


The elements i = 1,2,3 of the vector valued function 

f(c) are given by 


f^(c^) = B2 B^2l^sin(cr^ + (6°-6°)) - sin(6°-6°)] 

12(02) = \ ®3 Bi3Csin(o-2 + (6j-6°)) - sin(6°-6°)] 
and 

f^(0^) = B2 Ej B23[sin(c^ + (6°-6°)) - sin(6°-6° ) ] . (2. 51) 

The functional dependency among the state variables 
is apparent in their description. The dependency relations 
are given by 

X. - Xr- + 2^ = Constant 
X^^ - X2 + X^ = Constant 

These two relationships being functionally independent [65], 
the legitimate order of the state vector is 6-2 = 4, which 
agrees with our earlier observation that it should be 
(2n-2). These results can be easily generalised to the 


n-machine case 
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2.5.3 State model in the (n+m) dimensional state space. 
' (For non-uniform damping only) 


In this section a stole model in the ^^Cn+l) 
(=n+ni) dimensional state space for the n-machine case with 
non-unif ona damping is proposed. The state variables 
chosen for this purpose are 

, i = 1, 2, . . .n 

^i+n ~ ^ 0 ^ ~ ^ J ^ ^ ^ i— 1 j2,...iiij n,3 , 2 , . . . n 

(2.53) 


With these variables the state model is given by 


I = A X - B f (a) 

0 .= 0^ X (2.54) 


where 



"'^nn 

0 

nm 


B = 

J 

H 

I 

A = 


0 

o 

^ 0 



nm 

mm 



mm 



-j 





and 


0 = [0 


mn 


I ] 
mm-* , 


The f^( 0 j^)'s are the same as in equation (2.38). ¥e now 
exhibit overdescription in the above model by considering 
a 3-niachine example. 


3-iiiachine case^ 

As an illustration consider again a 3-iiiachine 
system with non-uniform damping. According to (2.53)? 
the state variable description is of the forms 
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=0)^ ; X4 = (63_-62) - (Sj 

X2=a,2 j I5 = (63_-63) - (Sj-e?) ( 2 . 55 ) 

X^ = J Xg = (62-^3) “ 








45 


The linea,r dependency among the variahles Z^ and Zg 
is evident. Therefore, the relationship "between them 
can he directly written e,s 

- Zg + Zg = 0 (2.57) 

Also no other functionally independent relationship among the 
state variables is found. Hence the legitimate order of the 
state vector Z is (6-1) = 5. This agrees with the earlier 
c onclusion ^that the minimal order is (2n-l) in the non- 
uniformly damped case. 

REMARK; Although the system models derived in Sec. 2.5.2 
and Sec. 2.5.3 are overdescrihed , they are, nevertheless, 
useful in obtaining a suitable decomposition of the system 
(2.1) into lower order subsystems as detailed in Chapter III. 

2.6 MIHIMAl ORDER REPRESENTATION OP AN n-MACHINE SYSTEM 
In this section the state models of an n— machine 
system with the state vector Z of the minimal order [lO] 
[(2n-2) for the uniform and (2n-l) for the non-uniform 
damping case] will be presented for the sake of completeness. 
These models are completely controllable and completely 
observable . 

Case (i) Uniform damping; 

The following state variables are chosen in this 
case with machine 1 as reference [10,69]. 



46 


= 0 ,^ - 03 ^ 



!1 

5 5 i = 2, 3, . = .n 

5 


With, these components the state model is 
X = X - f (a) 
a = X 


( 2 , 58 ) 


(2.59) 


where 


A.JJ - 


"■^^(n~l)(n-l) °(n-l)(n-l) 


^(n-l)(n-l) ^(n-l)(n-l) 


; B 


u 


t — T 

K / nM X 
n(n-l) nm 


0 


(n-l)m 


T 1 

^ t^®m(n-l) '^m(n-l)-^’ '^m(n-l) ^ 


^(n~l)(n-l) 

T 

"^(n-l) (m-n+l) 


and 


K. 


n ( n-1 ) 


-l(n-l) 


-I 


(n-1) (n-l) 


The row vector l 2 _(^ 1 ) Bas (n-1) elements each equal to 


unity. The fj^(cy^)’s are the same as in equation (2.38). 
The transfer function r 
the system is given by 


The transfer function matrix ¥-^( 3 ) of the linear part of 


= ' 


(2.60) 
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Case (ii) Ion-uniform dampings 

Here the state variables chosen are [69] 

X = (t) ; i =l, 25 ..,n 

^i+n " - ^i^ " " ^i^» ^ = 2,3,...n (2.61) 

With these variables the state model is 



2.7 COHOlUSIONSs 

This chapter has reviewed the state models of power 
systems in the light of some new concepts concerning un- 
controllability due to system overdescription. It is 
because of the fact that the state models in the 2n- 
dimensional state space both for uniform and non-uniform 
damping cases, the 2m -dimensional state space for uniform 
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damping and (m+n) state spa,ce for non-nniform damping 
case, exhibit redundancy in the state variable description 
either explicitly or otherwise, that uncontrollability 
in a mathematical sense is exhibited. This is explained 
in this chapter on the basis of physical dependency 
relations among the state variables. Prom these the 
legitimo-te order of the state vector is obtained. 



CHAPTER III 


STABILITY AHALYSIS HSIHG VICTOR LYAPHITOV PUHCTIONS 
3,1 IHTRODUCTIOH 

In tHis chapter a ne\<r approach to the stability 
analysis of large scale power systems using vector Lyapunov 
functions is introduced. Stability oanalysis of power 
systems using Lyapunov’s method has been carried out by 
several research workers [3~34] and is by now well 
documented in the literature. It is well known that the 
method obtains clearing times of faults in a one step 
integration of the swing equations thereby reducing the 
computer time considerably. Although encouraging results 
have generally been found [ 12513 , 17 » 20] , the modelling 
techniques had to be necessarily simple because of the 
difficulty in the construction of suitable Lyapunov functions 
with detailed models. As yet, the construction of Lyapunov 
function for multimachine power systems with models other 
than 'the constant voltage behind reactance' model is 
largely an unsolved problem. It therefore calls for new 
theoretic approaches to the problem. Two such promising 
approaches are 

(i) Dynamic equivalencing [53-64] 

(ii) Use of vector Lyapunov functions [39, 5l] 

The first approach consists in obtaining a simplified 
model which is then used to construct a Lyapunov function. 
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In the second method the composite system is decomposed 
into subsystems. Scalar Lyapunov functions are constructed 
for these subsystems. These Lyapunov functions together 
with the interconnections are then used to derive conditions 
for the stability of the composite system. The present 
chapter is devoted to this method of studying the stability 
properties of large scale power systems. 

The concept of vector Lyapunov functions was first 
proposed by Bellman [40] and Bailey [39] demonstrated its 
usefulness for studying the stability of a complex composite 
system using the idea of decomposition and the comparison 
principle [72]. The technique was restricted to a particular 
class of interconnections, i.e, linear interconnections only. 
However, it opened a new line of research in applying 
Lyapunov's method for lai^ge scale systems. Ilonlinear inter- 
connections were proposed by Piontkovskii et al [49l but the 
conditions were overly restrictive. Recent works' 'of 
Siljak [42-45], G-'amjic and Siljak [47,48], Thompson [46], 
Araki and Hondo [4l] and Miechel [50,5l] give' a good 
exposition of these concepts and enables one to apply the 
results to practical systems such as power systems. The 
work of Siljak [43] exhibits the nature of interconnections 
in a simple manner which nay be useful in -assessing the 
'connective stability* properties of complex systems. 
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This new technique of stability investigation of 
the complex system through subsystems is a two level 
concept. At the outset the composite system is represen- 
ted in the form of simpler subsystems and their interconnec 
tions. At the lower level, each of the subsystems is 
tested for its stability. A scals^r Lyapunov function 
satisfying the necessary sign definite properties is 
constructed for each of them in the usual manner. On the 
higher hierarchial level, these scalar functions are taken 
as components of a vector Lyapunov function [39] defined 
for the composite system. Using the constraints on the 
interconnections between the subsystems a vector matrix 
differential inequality is obtained. Applying the 
comparison principle [72] stability of the overall system 
is then ascertained. 


Applications of the vector Lyapunov function 
concept in physical systems in general and power systems in 
particular do not seem to have appeared in the literature 
so far. The power system, being a large complex system, 
is a natural candidate for using such a technique. This 
was hinted at in reference[35 ] .The main impediment in the 
successful application of the theory appears to be the 
difficulty in achieving an effective decomposition procedure 
of the complex system such that the subsystems have the 


desirable properties required for the construction of a 
vector lyapunov f inaction. Most 

. „ & 46SS3. 


Afic. Ne. 
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the' literature presuppose the existence of subsystems 
together with their interactions. But in a physical 
system an exactly reverse situation exists s a complex 
system is given and the first step is to obtain a suitable 
decomposition before applying the theory. 

This chapter therefore addresses itself to the 
problem of proposing a decomposition technique for large 
scale power systems and demonstrating the application of 
the vector lyapunov function approach for the stability 
analysis of the composite system corresponding to the 
post-fault state. A scalar lyapunov function of the lure'- 
Popov form is then obtained for the composite system from 
the vector Lyapunov function to estimate the stability 
region. 

3.2 SOI® ASPECTS OP VECTOR LYAPUl'TOV PIPICTI01TS[42-45,47,48] : 

Consider a continuous autonomous dynamical system ttS 
described by the vector differential equation 

i = f (2) (3.1) 

where Z is an n-dimensional state vector and the vector 
function f (X) satisfies the Lipschitz conditions so that 
the solutions of (3.1) exist and are unique and continuous 
for all initial conditions X^ of the state vector Z. 

Also f (0) =0 and the origin X = 0 is the unique equili- 
briiim point of the system (3.1). 
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Let the system 5 te decomposed into s dynairiical 
mb systems^, described by 


s 


Si =£i(Sl) + SijfiPf. i=l,2,---s 


j=l 




(3.2) 


where is the state vector of the subsystem of order n^, 

— i^— is a vector valued function of order n^, 
e. . are interconnection numbers to be defined later, 

_L. J 

_gij(y_^.) is a vector valued function of order n^^, 
is the output of the subsystem of order m^^,- 
and is a vector function of order 


It is assumed that the interacting function g. -(y^) and 

1 j J 

the function h^-(Z. ) are such that 

J J 


II |[ <ni. II X. 


(3.3) 


where r| . . are non-negative numbers and ||Z. || denotes the 

1 J J 

Eucledean norm of Z- given by (Xt X. 

o J u 

The state vector of the system S is given by 

S = [2i» ^2, ... Xg]^. (3.4) 

Accordingly the dimensions of and5.; , i=l,2,...s are 


related by 


n 


n 

.1 ^i • 
1=1 


(3.5) 


The subsystem (3.2) will be defined as ’forced subsystem' 
with the vector functions g.= .r(y.;j) representing the action 
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of the jth subsystem on the ith subsystem* If the summation 
terms on the right hand side of ( 3 . 2 ) are equated to zero, 
the subsystem description reduces to 


4 = fiUp 


Xl = h. (x^) . 


(3.6) 


'Ihen ( 3 . 6 ) defines the ' free' subsystem i=l, 2,...s. 

The interconnection numbers e. . assume values of either 

1 tJ . 

zero or one depending on whether or not -S acts on '5'.. 

J -L, 

Thus 


e^^ =1 if S j acts on^S^ 

e . . = 0 if 'S. does not act onS., 
J J 


(3.7) 


Using these numbers, it is possible to exhibit the inter- 
connections existing between the variotis subsystems through 
an ’interconnection matrix’,'^, and also illustrate the 
structural changes that can take place.. For example 
consider the system shown in Figiure 3.1. It consists of two 
subsystems a,nd of order, n^^ and n 2 respectively. With 
all the switches closed, the interconnection matrix is 
given by 


^ = 


1 

1 


1 

1 


indicating that all the interactions i,D = 1,2 are 

0 0 

on. Other forms of are possible such as |^q 0 
etc., depending on whether the corresponding g^^’s act or not. 

«L J 


o' 


0 1 

0^ 

9 

1 0_ 
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It will be ass-imed that . each of the free subsystems 
(3.6) are asymptotically stable. lor such’ systems, the 
following theorem holds [73]. 

THEOREM 3.1» let the eq.uilibrium of the system of 
differential equa.tions 

ii = (3.8) 


be asymptotically stable. Then there exists a positive 
definite decrescent function with a negative definite 

derivative. 


Proof of this theorem is contained in reference [75]. 
As a result of this theorem, satisfies the following 
conditions 


V. (xp < - II II ) 

where = (Grad 7^)^ is the total time 

derivative of 7j^(X^) along the solutions of (3.8), and 
Grad 7^ is the gradient vector given by 


(Grad 7^) 


I r Ifi 

Law * 'Sr . ^ 
^±2 


3X^ 




The Z- -‘s, 3=1, 2,,. n. are the elements of Z. , The 

X J X X 

functions ^ ±2 ^3 comparison functions 

belonging to class k . Some definitions in connection 
with theorem (3.1) a3:*e given in Appendix A for completeness 
and clarity. 
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In order to infer stability of the composite system 

it is necessary to establish conditions that the interacting 

functions have to satisfy. It will be therefore assumed tha.t 

there exist bounded functions C . -(S.) such that 

^ D t3 

(grad vp® fi.pxp < E, .(zp ®j3(l|Xj|| ) , 

i=l,2,...s. (5.10) 

Existence of these conditions crucially depends upon the 
nature of the Lyapunov function and the interacting 

fu.nctions. These in turn depend upon the way in which the 
system -S' is decomposed. 

The properties of given by (5.9) and the 

inequalities governing the interconnecting functions 
K- i(Si) given by (3.10) are now unified by taking the total 
time derivative V'j_(Zj^) of ’7^(Xj_) along the solutions of the 
'forced’ subsystem (^j^) . This gives 

7i(ZiI = ^^(2^) + (grad V^)^ gj, j(5j-) 


s 



i“lj2y...s . (3.11) 

Now define a vector Lyapunov function as 

and a comparison vector function ¥ as 

W = «23 ... 5^3]* . (?.13) 
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Using ( 3 . 12 ) and (3.13) in equations (3.11) the following 
vector matrix differential inequality is obtained. 

< R U ( 3.14) 

where Y^, ... Y^'f . (3.15) 


The elements r. • of the matrix R are given by 

1 J 

r - . =-6 j ^ ■ 




ID 


( 3 . 16 ) 


where 6^^ • is the kroncker delta symbol with 
6*1 . t= 1 for i = 3 

i J 

and 6t. = 0 for i ^ ji 

-J- J 

The are given by [48] 

-L J 

J. . = 6*. Sup ?. 4 (X.) + (l-^-J max(0,Sup ^..(X-) ). (3,17) 

Ij 1 J IJ J -I* J J J 


Finally to conclude stability of the composite system, the 
matrix R of (3.14) with its elements r^^ given by (3.16) ±i 
tested in terms of the following theorem [42,45,48]. 


THEOREM 3.2: The eouilibrium state X = 0 of the composite 
system S is asymptotically stable in the large if the sxs 
matrix R = (r- •) defined by (3.16) satisfies the conditions 

1 J 



^11 

^12 . 


... i-if 

‘^11 

ri 2 .. 

"■is 

H 

H 

A 

0 

^21 

^22 


^21 

• • 

"*22 ** 

• ♦ ♦ 

"*28 

♦ * 






■j 

"■si 

rg 2 .. 

r 

S3 


(3.18) 
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Proof of this theorem is contained in ref erence [48] , and 
is included in Appendix B . 

5.3 SOME PROPERTIES OE THE R-MTRIZ; 

1. Every element r- . of R is related to a single 

-A- J 

interconnection between two subsystems S. and 5'-. 

O 

2. The matrix R has non-negative off-diagonal 
elements. For such matrices the following lemma holds [74]. 

EEMMA 2,1; A real sxs matrix R (= with ,• >. 0 (for 

X J X J 

all i ,3 =l,2,..,s, ip^j) has all eigenvalues with 
negative real parts, i.e. 

Re(X*jj.) <0 , for all k = 1,2, .s, 

if and only if the inequalities (3.18) is satisfied. 

According to lemma 2,1, the matrix R satisfying 
(3.18), is Hurwitz. 

3. The matrix R is also negative definite. 

4. The negative of matrix R, i.e,, -R, has all 
off-diagonal elements non— positive , Such matrices are 
called ’Metzler' matrices [48,75] and have extensively 
been used in economic theory [76], 

3.4 SCALAR lYAPIMOV FIMCTIOI FROM VECTOR lYAPOTOV FDHCTIOl; 

In this section a scalar Lyapunov function Vg(Z) 
will be constructed for the composite system using the 
vector Lyapunov function. This scalar function is useful 
in estimating the stability domain for the overall system. . 



60 


Such, a region exists ■whenever the Lyapmiov f-unctions 
^i^— i^ of the subsystems do not satisfy conditions (3.9) 
in their entire state spaces but satisfy only in a 
certain region arouind the origins of their respective 
state spaces . 

¥e choose a f'unction V (Z) of the type 

( 3 . 20 ) 

as a candidate for the scalar Lyapunov function for the 

m 

composite system where d is a constant s-Tector (d^,d 2 ,..dg) 

with its elements 0. Sin'e every element of d is 

positive and every element of Y is positive definite, Y 

is also positive definite [48], Also, taking the time 

derivative of Y (X) yields 

s 

^ R (3.21) 

Since R has non-negative off-diagonal elements and 
satisfies (3.18) £ is Hurwitz. Hence for any given 
s-vector b ^ 0, it is possible to find the vector d of 
the form shoTinn such that 

b-^ = -d-^ R . (3.22) 

Therefore (3,21) gives 

VgCZ) < -b^ w (3.23) 

■vfhich shows Y(X) < 0. Therefore (3.20) indeed 
constitutes a Lyapunov function for the overall system. 
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When the estimates on the subsystem Lyapunov 
function i = l,2,,..s are linear (i.e. 

II -i II ^ ^ If 2,3 ;± = 1,2, ..s), it is possible 
to use the following theorem .[48] to obtain a scalar Lyapunov 

m 

function 'Vg(X) of the form H for the composite system. 


THEOREM 3.3* ^be ec^uilibrium state X = 0 of the composite 

system , S is asymptotically stable in the large if there 

exists a symmetric positive definite sxs matrix Q such 

that the symmetric positive definite sxs matrix H (= h. .) 

_L J 

which is a solution of the Lyapunov matrix equation 

R^H + H R = -Q (3.24) 


has all its elements specified as 


> 0 , 1^3 

>0 , 1 j . 


(3.25) 


Proof of this theorem is contained in reference [48]. 


Since the estimates on 7^ and 7^ are linear, it is 
possible to obtain a relationship of the form 

- ^i^ 15 » i = 1,2, ...s (3.26) 

where Ej[_'s are constants greater than zero. Then ¥ can be 
related to by the inequality 

7^(X)>K*¥ (3.27) 

where K* = Min i = l,2,.,.s, How choose a vector 


function Z 
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= H (3,28) 

where H is a positive definite symmetric matrix and 
consider a fimction of the form 

s "^i 

Vg(X) = 2 I /^(HV^)d7^ (3.29) 

i=l o *■ 

as a Lyapimov fimction of the composite sirstem. Expanding 


(3.29) one obtains 

V. 


V^CX) = 2[ / + ... + V3(x^)]d73_ 


0 

7. 


^ +li22’^2^^2^ + ••• + 1^28^8 ^^s^ 


7 


+ 


[h3i7,(Xp 


>^S2''2 


(X2)+...+ll^3Y3(X3)]aTg]. 


On performing the integrations, 


VgCX) = + 2h^2Vj_V2 + ... + hggYl 


(3.30) 


Olearly 7 (X) is positive definite, faking the time 
derivative of (3.30) 


VgCX) = 2 z’^dp k 

m 

= 2 t EY 

V ^ 

Using (3.14) we obtain the inequality 
< 2 7^^ H E ¥ 


(3.31) 
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and. in view of (3.27), 

YjX) < H R ¥ 

< K*[¥^R^ H ¥ + H R ¥] 

:< -K* Q I . (3.32) 

Since R is Hurvritz and H is positive definite symmetric, 
ty assumption, a positive definite Q always exists which is 
a solution of the Lyapunov matrix equation 

R^ H + H R = -Q . (3.33) 

Hence V (X) < 0. It may be noted that, the restriction 

on the elements h. . of H have been relaxed since the 
estimates on V^(X) and Y^(X) are linear. 

Since YjX) > O' and V„(X) < 0, V„(Z) is indeed a 
Lyapunov f-unction for the dynamical system 

¥ith this background, it is now possible to analyze 
the stability problem of the large scale power system. 


3.5 THE POWER SYSTEM PROBLEM; 

3.5.1 Decomposition; 

The swing equations characterizing an n-machine 
power system as indicated in equations (2.9) of Chapter II 
are given by 

d^6. d6- n 

^ rt- = 

i=l,2,...n . (3.34) 
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As stiown in Chapter II, Section 2.3, these equations can 
he cast into a vector matrix differential equation of the 
form 

1 = A X - B f (a) 

m (3.35) 

2 = C X 

where the state vector X and the vector o, and the 
matrices A,B and 0 are all of proper order depending on 
the state space in which the system is represented. 

Various forms of these models have been explained in 
Chapter II. A key to the application of the vector 
Lyapunov function is the ability to decompose the composite 
system into subsystems. The system models in the minimal 
state space i.e., (2n-l) for non-uniform damping and (2n-2) 
for uniform damping are not amenable to a decomposition 
where the interactions are functions of the state variables 
of the interacting subsystems as required in equation 
(3.2). Hence decomposition in the minimal state space 
was not possible. Ho v? ever the system descriptions 
discussed in Chapter II in the n(n-l) and (ii-tm) state 
spaces for uniform and non-uniform damping cases respec- 
tively lend themselves to an effective decomposition 
procedure. Therefore with these models, the decomposition 
procedure will be illustrated for a three machine system 
initially and will be later generalised to an n-machine 


case. 
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Case (i) 3-Eiachine system with, uniform damping 
( = X 5 i=l,2,3) • 

In this case the state model .as derived in 
Section 2.5.2 of Chapter II is of the form 



(3.36) 


with the state variables defined as 
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II 

— (O2 f 

-4 = ( V ^ 2 ) - 

( d °~ 6 p 


-r 0) 

■^"2 - 1 

3 

1 

= (6^^6^) - 

(6°-6p 


= li )^ 

-"3 i 

^6 = (^ 2 *^ 5 ^ - 

(62-6^) . 



The i=l,2,3 are the same as those defined in 

equation (2.59') of Chapter II. Decomposition into the 
3 subsystems is achieved by choosing the state variables 
from equation (3.57) as follows s 


SUBSYSTEM 1: Choose 

X-, =0)^ - uj 


(5.38) 




as the state variables for tiiis subsystem. Extracting 
equations corresponding to those states from (3.56), we 
get 




\ 

i 


—X 

0 

■i 

^1 


’1 + 1 ’ 
1^2 

• 





1 

0 

i 

.^4_ 


0 

u J 


° J 


(T-, = [0 



1 

M 2 


_o 

1 ] 1 


i 

-^ 4 . 


f ^(a^) 


(5.59) 
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SUBSYSTEM 2: Choosing 


^2 = ^2. “ “•= 


I, = (62 - 6^) - (6 


6 °) 


(3.40) 


as the state variables of this subsystem, and extraction 
of the state equations corresponding to these variables 
from (3.36) results in the description of the forced 
subsystem 2 as 


^2 

• 


-X 

0 


X2 


4 + i 

«i % 

Q 

ro 

1 

1 

J 



1 

0’ 


X 5 


0 


0 

- - 


q(cri) 


1 

M, 

0 


^2 = [0 


f ^(0^) 


X. 


(3.41) 


SUBSYSTEM 3; Adopting the same procedure as for the first 
two subsystems the state space description of this subsystem 
with the state variables 



Se = (^2 " ^3^ ” ^^2 " ^3^ 


(3.42) 


will be of the form 
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^3 


-X 0 

rr 

Xx 

3 


■ 1,1 


1 

“2 

0 

«=; 







1 

M 

CT\ 

1 


1 0 

1 


0 

L- J 


0 

L 




1 

0 


£3(^2) 


(3.43) 


= [0 1 ] 


X, 


The system (3.36) is thus decomposed into 3 subsystems 
described by (3.39) » (3.41) and (3.43). These ecxuations 
together constitute the composite system analogous to (3.2) 
and thus represent the dynamics of the 'forced' subsystems 
1, 2 and 3 respectively. The last two terms on the right 
hand side of the equations (3.39), (3.41) and (3.43) are 
identified as the interactions among the subsystems. The 
description of the 'free' subsystems will then be as 
follows; 


SUBSYSTEM 1: 


' H 

• W 

i 


-X 

0 


1 


ri + 

:m/ 

1 -l 

1 

• N 

i 


' 1 

0 


Y4J 


0 



(Tl = [0 1] 


(3.44) 
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SUBSYSTEM 2l 


'^2 


X 

o" 


^2 


-1,1-1 
M3_+ M^ 


j 

1 

0 

i 

i 



0 

L 


and 




(5.45) 


SUBSYSTEM 5: 



- 

-X 

0 

1 

1 

E3 


rl . 1 -1 
M2'^ 



^1 

0 

j 

.^6_ 


f 

0 


a 


3 


[0 



'~<r 


6 


(3.46) 


It is evident from the above equations that the free sub- 
systems represent the dynamics of a hypothetical 2-machine 
system cast in the lure ’-Popov form. Thus subsystem 1 
consists of ma,chines 1 and 2, subsystem 2 includes machines 
1 and 3 and subsystem 3 comprises of machines 2 and 3. 

It may be noted that the decomposition is a mathematical 
one and not a physical one as in Eron's 'tearing technique’. 
The interactions among the subsystems take place thro^ugh 
the second and third terms on the right hand sides of 
equations (3.39), (3.41.) and (3.43). These are nonlinear 
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due to nonlinear nature of occurring in them. The 

interconnection numbers e . . exhibiting these interactions 
among the 3 subsystems are given by the following inter- 
connection matrix 


(1) (2) (3) 


(1) 011 

(2) 1 0 1 (3.47) 

( 3 ) |_1 1 0 

This matrix is invariant in structure. A block diagram 
representation of the interconnected system S, with the 
subsystems and their interactions explicitly displayed 
is shown in Figure (3.2). 


Case (ii) 3-machine with non-uniform damping 
(^j_7^^j>i»j=l>2,3)" 

The state model of a 3-ma.chine system obtained 

in Section 2.5.3 of Chapter II in the (m+n) dimensional 

state space with the state variables defined by 


11 

B 

H 

X4 = (6^ 

- 63) 

0 H 
<0 

1 

- 

X2 ^ “ 2 ; 

X5 = (61 

- '63) 

- (6° 

- 6°) 

0 

3 

II 

M 

Xg = (62 

- 63) 

- (6° 

- 


(3.48) 


is of the form 




Subsystem 3 

FIG.3.2. SUBSYSTEM INTERCONNECTIONS IN A 3MACHINE 
POWER SYSTEM (UNIFORM DAMPING) 
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“^1 

0 

0 

0 

0 

0 ' 



^2 


0 

-^2 

0 

0 

0 

0 


X2 


= 

0 

0 

-X3 

0 

0 

0 


^3 

• 


1 

-1 

0 

0 

0 

0 


♦ 

m 

1 

1 

0 

-1 

0 

0 

0 

i 

« 

* } 

0 

1 

-1 

0 

0 

0 

• j 

_"^6_ 




0 


0 

0 


1 

1 



1 1 

D D 

0 0 

0 0 

0 0 


f 2 ^ ^2 ^ 


(3.49) 



02 



0 0 0 
0 0 0 
0 0 0 


1 

0 

0 


0 0 

1 ; 0 

0 1 


X. 


1 


."6. 


Direct decomposition of (3.49) into subsystems described 
by equations of tiie type (3.2) where the interacting 
functions g. .(y,. .) are functions of the state of the jth 
subsystem only is not possible. This difficulty is overcome 
by augjt. 3 nting the first three state equations 

^1 "^1 ^1 *■ 
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♦ 1 

^2 ~ ~^2 ^2 




Effectively, this implies that the state description of 
the composite system is in the (&-3y^ dimensional state 
space. This does not, however, affect the results and is 
only an artifice to achieve the desired decomposition. 

Decomposition into three subsystems is now 
stra-ight forward and is carried out in the following way; 


SUBSYSTEM 1: The state variables selected for this 
subsystem are 


Z2 = <^2 (3.50) 


Extracting the corresponding equations from (3.49), 
subsystem 1 is described by 






0 




^ 3 ( *^3 ) 



X.. 


(3.51) 


a = [0 0 1] 
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SDBSXSDEM 2: Here the state variables are 

Xi = 

= ( 0 ^ ( 3 , 52 ) 

X5 = (61 - 63) - (6j - 6°) 
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CT^ = [0 0 1] 



■with, the state ■variables given by 

Xg = (62 - 6^) - (6° - 6^) 


(3.54) 


(3.55) 


Thus the 3-niachine system is decomposed into 3 subsystems 
described by (3.51), (3.53)and (3.54). Each of these 
represents a 'forced’ subsystem analogous to the system 
( 3 . 2 ), Again with the last two terms on the right hand 
sides of these equations identified as 'interactions', 
the 'free' subsystems are described by; 


SUBSTSTM 1; 




t 


1 T 

“1 : 


^2 

- 

1 

-Mj 


-^4- 




0 

{ 




[0 0 1 ] 



(3.56) 


SUBSYSTEM 2: 


'h 



0 

0 ' 


Zi 


-1 - 
^1 


V 



0 

-X 

3 

0 


z, 

3 

— 

1 

-M 3 

^ 2 ^ ^ 2 ^ » cr 2=[0 0 l ] 

Z 3 



1 

-1 

0 


Ys- 


0 _ 


1^5 J 


(3.57) 
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and 

SUBSYSTEM 3^ 




-X 2 

0 ' 

0 


X 2 


1 

M 2 


r~ 

^2 


=Z 

0 , 

3 

0 


^3^ 

- 

1 

” 

3 

; 0 ^ = [0 0 1 ] 


^6 


1 

-1 

0’ 


^6 


0 

i 



(3.58) 


Bach, of these free subsystem models can be identified with 
that of an equivalent 2 machine system with non-\miform 
damping. The first comprises machines 1 and 2, the second 
1 and 3 and the third 2 and 3. A block diagram represen- 
tation of the composite system exhibiting the subsystems 
and their interactions is given in Eigure 3.3. Einallys the 
interactions represented by the last two terms on the right 
hand sides of (3.51), (3.53) and (3.54) are exhibited by 
the interconnection matrix CDof equation (3.47). 

The procedure of decomposit in described will now' 

be generalized to an n-machine system for both the uniform 
and the non-uniform damping cases. 

Case (iii) n-machine system with uniform damping; 

State model in the 2m-dimensional state space for 
the n-machine case was derived in Section 2,5.2 of Chapter 
II. Decomposition of this system into 'm' subsystems is 
carried out in a manner outlined for the 3-machine case. 
Then a typical ith subsystem involving machine 'k’ and 'o’ 
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is of the fora 


n 


X. = A. Z. - B. f.(cr.) + T e. B f (cr ) 
-1 1 -1 ”1 i'' iq -q ~q^ q-^ 




a. = 4 


i=l, 2, . . ,n 


( 3 . 59 ) 


where X^ is the subsystem state vector. The matrix Aiand 

vectors C., B. and B are obtained as 

—1 ' —1 — q 


^1 = 

-X 

0 


j. 

Cj_ = [0 1]^ ^ = ; 

i + i 


1 

0 ; 


1 

' 0 ^ 



1 1 


1 

and B = + 

. \ 

or 

®3 

-q 



; 0 ' 


^0 , 


( 3 . 60 ) 


As before the numbers e. constitute the elements of the 
interconnection matrix"^ . This matrix is of the form 

= (K^ E - 21 ) , (3.61) 

' v^-m vrm mm' ” 


mm 


nm nm 


mm 


where I is an identity of order m. 

The free subsystem 6"^ is finally described by 


T 


Xq = X. - B. f . (o.) J a. = 0 . X. 


( 3 . 62 ) 


and represents the dynamics of a 2-machine system, 

Oase (iv) n-machine system with non-unifoim damping: 

The state model of the n— machine case with non— unif oim 
damping in the (m+n) dimensional state space developed in 
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Section 2.5.3 will be used here to obtain a suitable 
decomposition into ’m’ subsystems. Adopting a procedure 
ideiitica.1 to the 3—machine case, a typical ith subsystem 
involving the kth and jth machines (lc< 3) will be as 
follows. 


^i 


A. X. - B. f. (a. ) + y e, B f (a ) 
1 -1 -1 iq -q q^ q^ 


(3,63) 


- T 

Si 


a 1 , 2 , . . . m 


where X^^ is the subsystem state vector of order 3. 5)he 

matrix A. and vectors G . , B. and B are: 

1 —1’ —1 — q 


A, = 



0 

(T 


r 1 

Mk 

0 

-'3 

0 

: 1 B . = ; 

1 

■■“3 


-1 

0 ^ 


0 . 


Si = [ 0 


if 


and B = + 

-q 


1 


f*- 1 

0 

1 

% 


0 

j 

or + 

M. 


L ’ 

J . 

,1 0 


0 ( 





(3.64) 


The interconnection structure among the subsystems is 
exhibited by the interconnection matrix '^given by (3.61). 
lastly the free subsystem is of the form 


Si = h Si - £1 fi(m) 


S Si 


( 3 . 65 ) 


which again represents the dynamics of a two machine system 
with non-uniform damping. This completes the decomposition 
procedure. 
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3.5.2 Construction of Lyapunov function for subsystems; 

For tile free subsystems described in Section 3.5.1 

scalar Lyapunov functions will noxir be constructed. Eaeb 

of these subsystems is in the Lure '-Popov form with a 

single nonlinearity. The pairs (A^, B^) and (A^, 0^) are 

easily verified to. be completely controllable and completely 

observable. Hence the triplet (A^, C^) constitutes a 

minimal realization [69] for the hypothetical 2-machine 

subsystem. The transfer function of the linear part of 

these subsystems are given in equations (2.32) and (2.28) 

of Chapter II both for uniform and non-uniform damping 

cases. Techniques for generating Lyapunov functions for 

such systems are well known [7~10,14,15] . Here we use the 

Moore-Anderson theorem for constraoting the scalar Lyapunov 

function for each of the subsystems. This theorem is 

stated without proof in Appendix C. Following the method 

of reference [lO], a Lyapunov function F^(X^) of the type 

^i 

TTlUP = ^ Ti + 2Pi I fi(spas. (3.56) 

O 

is constructed in the respective state spaces of the 
subsystems <5^, i=l,2,,..s. 
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with. |3 . = ~ . 

1 X 

the matrix 


Similarly for the non-unif o rm damping case 
takes the form 


P.= 

1 




r"k 


0 




i. 


xyv5-r 


0 




% “l 


SlIi 


( 3 . 68 ) 


and p. = (~ + - ) . 

^ '^k '' j 

The time derivative V(X) of this lyapunov function is 
given by [10 ] 


Tpxp = - - 2fpap <.. (3.69) 

and negative definite for all lying in the open interval 

rn 

n-af j_) and the positive semidefinite matrix 
is obtained from the solution of the Lyapunov matrix 
equation 

A^P+P1=-LL^. (3.70) 


This implies that the subsystem is asymptotically stable 
in the region of interest. It is precisely this property 
that has resulted in a successful application of the vector 
Lyapunov function to the stability analysis of the composite 
power system. 

3.5.3 Stability Regions: 

The utility of Lyapunov functions in power 
system stability analysis lies in the computation of 
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their stability domains. In tMs section methods of 
determining these regions are discussed. 


The nonlinearities occurring in power system models 
are known to satisfy the Popov’s sector conditions 

0 < crj_ 1 “ » i = 1^2, . . .m 

/ , (3.71) 

fiCCi) = 0 for = 0 

for only those that lie in the closed interval 


-bi 

< 


where 

= tc+2(6® - 

(3.72) 

’ai 

= x-2(62 - . 



Therefore a region of stability exists aroimd the origin 
of the state space of the subsystem Hence for the 

composite system too there exists a region of stability. 

This region can be estimated by a knowledge of the stability 
regions of the individual subsystems [77]. Methods are 
available to compute these regions [67,68], Alternatively, 
the scalar Lyapunov f motion constructed in Section 3.4 
can be used to estimate the stability domain for the 
composite system directly. While the former technique 
needs further investigation, the latter procedure is 
adopted here to compute the stability regions of the 
composite system. A recently proposed method [3l] is 
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applied to estimate the nearest unstable equilibrium 
point for this purpose. 

• We shall now illustrate the complete procedure 
with reference to a 3-iaachine system. 

3.6 APPLICATION TO A 3-mCHINS SYSTEM: 

3.6.1 Decomposition and application of vector 
lyapimov method: 

The example chosen to illustrate the applica- 
tion of the vector Lyapunov functions is the 3-niachine 
system (Figure 3»‘^i-) of reference [5]. The system data 
is given in Table 3.1, and the transfer admittances 
between the internal nodes of the machines are as 
f ollows : 

^12 = 0 + p.u. 

^13 = 0 + (3.73) 

and y2^ = 0 + j3 p.u, 

A disturbance is simulated by changing y^^ from 3 p.u, 
to 1 p.u, and assigning arbitrary velocities to machines. 


Table 3.1 - SYSTEM DATA 


H 

Voltage 
(p.u, ) 

Me ch. Power Input 
(P„i) Cp.u.) 

Inertia Constant 
(p.u. ) 

1 

1.0 

1.5 

0.02 

• 2 

1.0 

0.0 

0.002 

r 3 

' 1.0 

: -1.5 , 

; 0.03 ■ 
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For this composite system the post-fault stable 
eq.uilibriuni point is given by [ 5 ] 

6 ° - 6 ° = 18.78° 

6 ° - 6 ° = 58 . 86 ° (3.74) 

and 62 - 6 ^ = 40.08°. 

With \miform damping and X=1 a.ssumed, the state represen- 
tation of the system according to equation ( 3 . 36 ) is 




Decomposition yields the following 5 forced systems: 
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SUBSYSTEM 1: 


• 

-1 

0 ” 


^550 


11 

1 

0 

:si - 

0 ; 



L. 

- 1 - 


—3 



= [0 

iJSi 



SUBSYSTEM 

2s 




• 

-1 

0 


~83.: 

33 ' 

^2 = 

. 1 

0 ' 


0 



— 

— 

i 




50 

0 




500 
0 

(5.76) 




f 2 ( ^ 2 } 


50 

0 




33.33 

0 

(5. -77) 




02 = [0 1] X2 

and 

SUBSYSTEM 3: 

-1 0 
1 0 


% 


% - 


533.3 

0 


f^Ccr^) - 


500 

0 




33.33 

0 

(3.78) 


± 2 (^ 2 ) 


= [0 1] Z5 . 


The ’free' subsystems are given by 
SUBSYSTEM Is 

Pcmn 

! CTJ^ = [0 1] % (3.79) 


♦ 

-wl 0 i 


550 

s = 

: 

il - 


”1 

1 0 ■; 


; 0 : 


SUBSYSTEM 2: 

-1 0’ 


X 2 


1 0 


X 2 - 


83.331 
0 


J c72 = [0 1] X 2 (3.80) 
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and 

SUBSYSTEM 3: 


% = 


-1 

0 

X — 

533.3 

1 

0 

% 

0 


£3(^3) ; (J3 = [0 1] Z3 (3.81) 


In equations (3 •15) - (3.81), the nonlinearities 


f j_(cr^), i = 1,2,3 are 


f3_(cy^) = 2.0 [sin((j3_+(6°-6°)) - sin( 6 °- 62 )] 

± 2 (^ 2 ^ =1.0 [sin(<r2+(6j-6°)) - sin(6°-6p] (3.82) 

^^(o’^) =1.0 [sin((j^+(62~6^) ) - sin(62-S^)] . 


Tliese nonlinearities satisfy the Popov's sector condition 
for cr^, i = 1,2,3 lying in the closed interval 

-3.7974 < c73_ < 2.4858 

-5.1970 < a2 < 1.0862 (3.83) 

-4.5412 < <1.7420 . 

Also it is easily verified that the origins of the free 
subsystems (3.79) - (3.81) are asymptotically stable. 

Hence for each of these subsystems Lyapunov function of 
the type 

Vi(lj_) = Pj_ Xj_.+ 4 / fj_(s^)ds^, 1=1, 2,3 (3.84) 

o 

are con struct ed* The matrices for each of these functions 
are ^ 

Po. 00364 0.00182 I 


0.00182 0.00182 



88 


P 


2 


and P^ 
D 


0.024 

0.012' 

0.012 

0.012 

0.00374 

0.00187~ 

0.00187 

0.00187 


(3.85) 


We choose a function i = l»-2y3 and 

obtain the following linear estimates on [42, .78] 

0.0263 |{ X^11< v^(Xi) < 2.0 II X^j| 


0.068 II < V2(X2) < 1.41511 X 2 II (3.86) 

and 0.0265 II X^ II ^ ^^(X^) 1 1*415||X3|| . 


The derivatives v. (X. ) of the scalar functions v. (X. ), 

1 jL «i, -L 

i = 1,2,3 along the solutions of the free subsystems are 

^l(Xi) = - 2 ^ (0.00364 Xj + 2f^(o^)a^) = -$3_3( II • X^ll ) 

"" “ “i^^(0*024 X| + 2f2(a2) = ~«'23^ II - 2 II ^ 

^3(13) = - “17^(0.00374 X^ + = -§33(11 %|| ) . 

(3.87) 


The ■^3*3^ i = 1,2,3 will be identified as comparison 
functions. 

HOTE^ The derivatives ^j_(X^) of (3.87) are negative 
definite for cr^,i=l,2,5 lying in the open intervals 
(- 3 . 7974 , 2.54858), (- 5 . 1970 , 1.0862) and (-4.5412, 


1.7420 ) respectively. 
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Now the derivatives v^(Xj^) , i = 1,2,3 of v^CX^^) 
along the solutions of the forced subsystems are computed, 
This gives 

=-^h) +-^(2Xi+^4>*3^'"3> 

’<_-U 3 ( I SJ) + Uz "zs'll S2II ) n 3 % 3 ( 1123 II '> 


Similarly 

£ 2 ^- 2 ^ ^r^i IIX^II ) -$23^1^2 11^ ^23 ^3^n~3li ^ 

and 

^3(1^) ^ ^31 ^13^11 Sill ) + ^32 ^ 23 ^ 11^211 ) ^33^ II S3 II ) 

( 3 . 88 ) 


where the r. .’s are constants >_ 0, 
a 3 

Define the vector Lyapunov function 




(3.89) 


1 i. 

and the comparison vector function ¥ = [$i^,$23> ^^3^ 

Then according to equation (3.14) 

» 

v^< R ¥ where 
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Computation of 's in the above equation is a 

critical part in this procedure. Majorisation procedures 

as used in ffrujic and Siljak [48] gave rise to extremely 

conservative results even to the extent of indicating the 

system to be unstable although it is known to be steady 

state stable. It was threfore necessa,ry to consider the 

physical aspects of the problem in arriving at sharper 

estimates on the r. .'s. As is done in Gless [5], different 

^ J 

initial velocities were assigned to the machines with the 
pre-fault rotor angles as the initial conditions. For low 
velocities the system was found to be stable with the 
assumption that the velocities decrease monotonically , 

For larger initial velocities the system was, hovrever, 
unstable. ¥e thus compute the ©"^aluating the 

right hand sides of (3.88) at the initial point. Table 3.2 
summarizes the effect of the initio-1 values of the state 
variables of the individual subsystems on the stability 
of the composite system. The state of the system when 
different initial conditions 3.re assigned to state 
variables, is given in Column 5 and is compared with the 
state obtained in Table I of reference [5]. This compari- 
son exhibits a tendency towards instability by this method. 
Thus the system stability domain is decreased with the 
present technique. Further numerical experimentation is 
called for to increase this region. One of the reasons 
for such a conservative estimate of the domain seems to 



Table 5.2; Bffect of Initial Conditions on Stability 


S.No. 

Value of 
^1 

Value of 
^2 

. 

Value of 

j System state 

1 

1 

0.03 

0.09 

0.06 


Stable 


2 

0.00 

0.10 

0.10 


Stable 


3 

0.025 

0.075 

0.05 


Unstable 


4 

0.001 

0.0705 

0.0685 


Stable 


5 

0.02 

0.08 

0.06 


Stable 


6 

0.01 

0.07 

0.06 


Stable 


7 

0.01 

-0.05 

-0.04 


Unstable 


8 

0.04 

0.01 

-0.03 


Unstable 


9 

0.06 

0.05 

-0.01 


Unstable 


10 

0.01 

0.09 

0.08 


Stable 


11 

0.100 

0.00 

-0.10 


Unstable 


12 

0.03 

-0.02 

-0.05 


Unstable 
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and 



(3.93) 


•where K* is a constant greater than zero. The matrix 
Q as a solution of the Lyapuno"^ matrix equation 

R^H+HIl = -Q 


is given hy 


Q = 


1.00 - 0.256 - 0.932 

- 0.256 1.00 - 0.112 

- 0.932 - 0.112 1.00 


(3.94) 


and is positive definite and symmetric. Hence V_(X) is 
a lyap-unov function for the composite system. This 
function can be used to compute the stability region for 
the 3-niaohine system which is done in the next section. 


3.6.3 'Comparison with other Lyapunov f'unctions 
and stability regions: 

The Lyapunov function 'V_(X) in (3.92) will now 

be compared with the Lyapunov function constructed in the 

minimal state space. For a comparison to be valid both 

the functions have to be in the same state space [24]. 

To do so the Lyapunov function 'V„(Z) will be recast in 

s 

the minimal state space utilizing the linear dependency 
among the state variables in X. Considering the state 
variables X^^, X 2 , and X^ to be the linearly 
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independent components, the fmiction 'th.e minimal 

state space would be 




3 

Sm + \ / £ 4 ( 3 . )dS. 

-mi -m -m J 11 1 

1=1 o 


(5.95) 


where “ 1 “ ^ ^ 


T 


and 


P 


m 


0.001845 

-0.000935 

0.0009225 

-0.0004675 


-0.000935 

0.006935 

-0.0004675 

0.003475 


0.0009225 

-0.0004675 

0.0009225 

-0.0004675 


-0.0004675 

- 0.003475 

-0.0004675 

0.003475 


(3.96) 


Direct computation using the procedure of reference [lO], 
yields a Lyapunov fimction 7^(^) of the foim 
m 3 "^i 

^ 4 + . I / 

1=1 o 

with 



0.000952 

-0.000579 

0.000321 

-0.000193 


- 0.000579 

0.00635 

-0.000193 

0.00217 


0.000321 

-0.000193 

0.000321 

-O.OOOIOM- 


-0.000193 

0.00217 

-0.000193 

0 . 00217 . 


(3.98) 


and the ^^(^i)'^ are the same as in (3.91).’^* sbru.c'tu.te of 

is g(.vs.»s> €«. A|»|as.Tux£* C. 


Stability regions will now be obtained for the 


system using and '‘7^ by the technique of reference [3l]. 
The estimated unstable equilibrium point nearest to the 
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stable one according to this method is 

^ = [0 0 0 1.0862]^ (5.99) 

The values of the functions and Y^C^) computed 

at this point d.efine the stability regions. Thus 

0.591 (5.100) 

is a stability domain using (5.91) and 

\(^) < 0.588 ( 5 . 101 ) 

is the stability regions using (5.97). It cannot, however, 

be concluded that (5.100) yields a better result since 
the lyapunov fimctions differ in their quadratic part. 

From the entries in the matrices and P^^ it can be 
observed that tho trajectories reach the boundaries faster 
using the Lyapunov function '^^(X^) . It is therefore 
conjectured that vector Lyapunov functions method yields 
conservative results than the existing methods. But this 
disadvantage is outweighed by the other advantages of 
the technique. 

5.7 COKGLUSIONS; 

This chapter has demonstrated the applicability of 
the vector Lyapunov fionction to the stability analysis of 
large scale power systems. The key to a successful 
application of the method lies in an effective decomposition 
of the power system. In this chapter a decomposition 
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procedure is developed both, for 'uaiform and the non- 
uniform damping cases. This enables one to construct 
lure’ -type Lyapunov functions for the subsystems. Using 
these functions a scalar Lyapunov function for the 
overall system is obtained via a vector Lyapunov 
function. This scalar function is compared with the 
existing ones. Although it yields conservative results, 
the possibility of including more system details makes 
the new approach an attractive tool for a detailed 
stability analysis of the large scale power system. 



OHAPTBR IV 


TMMSIMT STABIIITY REGIONS FOR I'niLTIMAGHINE 
PO¥BR SYSTEMS 


4.1 INTRODUCTIONS 

One of the major impediments in the successfiiL 
application of lyapanov's method to practical systems 
is the computation - of the transient stability region 
around the post-fault stable equilibrium point. The 
estimation of this region so far has been carried out by 
computing the unstable equilibrium point closest to 
the post-fault stable one and evaluating the Lyapunov 
function V(X) at this point. This determination is 
based on, the following geometrical reasoning. The surfac 
V(Z) = e where e is a positive constant surrounding the 
origin Z = 0 are closed. As the value of e increases 
these closed surfaces expand such that V(z) = is 
contained in V(X) = ^2 e2>‘e2,* However, when the 

sxirface V(Z) = constant passes through an equilibrium 
point these surfaces are no longer closed. At this point 

V(Z) = [Grad V(X)]^ X = 0 (4.1) 

since X is equal to zero at an equilibrium point. The 
region bounded by the closed surface 


7(X) = V(Z^) 


(4.2) 



98 


where is the ■unstable equilibri'um point closest to 
the origin, i.e. the stable equilibri'um point, is the 
region of stability. The determination of the equili^ 
bri-um point of the post-fa.'ult system invol'ves the 
solution of 

Il2(S) = 0 (4.3) 

■where is the set of equations on the right hand 

side of equation (1,3), These equations are nonlinear 
and have (2^-2) solutions. Of these (2^”^-!) are unstable 
solutions and one is the post-fault stable solution. 

While the determination of stable equilibrium is not 
difficult, the computation of the closest unstable 
equilibrium point is a formidable task, and involves 
considerable amcunt of computer time. Various minimi- 
zation procedure such as the steepest descent [6], 
Davidon-Pletcher-Powell [12], the steepest ascent [27] 
and Uewton-Haphson [28] have been employed for the 
determination of these unstable equilibrium points. Some 
reduction in computer time was achieved by luders [17] 
by approximating these points. Willems and Willems [16] 
defined a 2m -dimensional boundary based on the sector 

violation of the nonlinearities ard obtained a minimum of 
V(Z) on the surface. A similar technique was proposed 
by Zaslavskaya etal [29] and Ananda Mohan [22]. In all 
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these cases the problem involves the use of some sort 
of a minimization technique for obtaining a minimum of 
V(x) on the boundary. Bergen otal [34] and Di Oaprio 
etal [30] use series expansions and obtain a somewhat 
conservative bound on Y(x) . Some physical energy-type 
considerations were adopted by ¥illiams etal [20], In 
spite of these techniques, the computation of the region 
of stability has been one of the discouraging factors in 
the practical application of the Lya,punov theory. Recent 
worh of Prabhakara and Bl-Abiad [31] may be considered 
as a powerfifL contribution in this area and enables one 
to apply Lyapunov methods to practical systems. The 
method yields approximate values of the region of 
stability within reasonable error bounds. Using the 
analogy of a single machine connected to an infinite bus, 
the method approximates all the unsta.ble equilibrium 
points. In this chapter a theoretical basis to this 
technique is provided. The method utilizes the results 
obtained in connection with the determination of stability 
regions for multinonlinear systems as a generalization 
of the work due to Walker and Mc'Clamroch [68] and 
Weissenberger [67]. An altoamate method that gives 
improved stability regions based on a conjecture due to 
Murthy [24] is also described for practical implementation. 
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4.2 REG-IOFS OP ATTRACTION POR MULTINONLINBAR SYSTEMS i 
4.2.1 Prollen formulation; 

The problem of determining the regions of stability 
for systems cast in the Iirre’ -Popov fom consisting of a 
single nonlinearity was analyzed by ¥alker and Mc'Clanroch 
[68] and V/eissenberger [67] using the sector violation 
information of the nonlinearity. In this section the 
method will be extended to a mult inonlin ear system. 

Consider a system S described in an n-dimensional 

state space and of the form 

X = A Z - B f(o) (4.4a) 

o = C^ X (4.4b) 

where A, B and C are matrices of order nxn, nxm and 
nxm respectively, f C^) ^ vector valued function of 

order *m' similar to equation (2.38). (it may be noted 
that m ^ n(n-l)/2 as in tho power system problem). The 
ith component of f (o) is a function of tho iih component 
of the m-vector a, X is an n-dimensional state vector. 

It is assumed that (A, B, C) constitutes a minimal 
realization of the transfer function matrix ¥(s) of the 
system (4.4) where 

¥(s) = C^(sl - A)"^ B . (4.5) 

Further, let there exist constant matrices a«=diag(a|_) 
and Q=diag(qj^), c££> 0, 0, > 0, such that 
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tile function 

Z(s) = (a + Qs) ¥(s) + aG~^ (4.6) 

is positive real [ 8lJ guranteeing the existence of a 

lure '-typo Lyapunov function of the type 

^ a. 

m ^ i 

Y(X) = r P Z + J 2q^ j fi(s^)dSj_. (4.7) 

i =1 o 

valid in the sector 

5, i=l,2,...m (4.8) 

and ^i^^i^ ® ' ^i=® ^ 

In (4.6) the matrix G- = diag(gj^), i=l,2,...m. The 
nonlinoarities ^j_(h^) assumed to satisfy the Popov 

sector (4.8)for only those that lie in the closed 
interval 

cTj^— , i=l,2,.,.m (4.9) 

as shown in Pigure 4.1, k^ and £ ^ are constants greater 
than zero. letting 

2q^ = , 7(X) can he written as 

m m ^i 

V(Z) = Z-^ P Z + j; p. / f . (s. )ds. . (4.10) 

i=l ^ o 

Since the nonlinearities fj^(cT^), i = l,2,.,.m 
violate the Popov’s sector condition (4.8) outside the 
interval defined in (4.9), asymptotic stability of (4.4) 
in the large cannot be concluded. Therefore there exists 
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a region of stability around the origin of the state 
space. An estimate of this domain of attraction defined 
by the following inequality is to be obtained, 

i.e. , , V(X) < e (4.11) 


where e is a positive constant to be determined. Clearly 
on the hypersurface defined by 

V(X) = £ (4.12) 


and its interior, the Lyapunov function and its time 
derivative 'V(X) satisfy the necessary sign definite 
properties namely, 


■v(x) •> 0 

V(X) < 0 

together with grad V(X) ^ 0 
where grad V(X) = [3/9X^, 9/3X2 


(4.15) 


9/9X^]^ V(X) 


4.2.2 Determination of the region of attraction; 

The .procedure to obtain an estimate of the region 
of attraction using the sector violation information of 
the nonlinearities will now be discussed. 

Let the matrix C be partitioned in the form 



^11 

C12 

^15 , * * 


C -” [ f 0^2 y • • ^ ^ ~ 

O21 

• « 

^22 

• • * 

^23 ** 

« • « • 

°2m 

• • 


Cnl 

C 0 
n2 

^n3 * • 

C 

nm 


(4.14) 



Then ec^uations (4.4'b) and (4.14) yield 
T 

= C X , 1=1, 2,... m. 

Row define 2m hyperplanes in the following ways 
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(4.15) 



Z =k3_ , 

i 

““ 1 ^ 2 ^ • # • m 

^3 



~ 1, ^ 2 ^ « o m 

in the n-dimensional 

state space 


(4.16) 


T 

bonndary of the closed region { Xs 0^ X < 

i = l,2,...m}. For example, in a single machine connected 
to an infinite bus the boundary JB in the 2-dimensional 
state space is a set of parallel lines. OnIB and its 
interior the 7-function and its time derivative satisfy 
the necessary sign definite properties. Hence the 
region contained by IB is in the exact stability domain. 

A search for a minimum of 7(X) on this hypersurface, 
therefore, gives an estimate of the region of attraction. 
We thus find the largest of the V-surfaces that still 
lies inside the surface IB . 

In order to establish this region consider a 
typical ith hyperplane of (4.16) given by 


<^1 = ^ 5 = t . 


(4.17) 


In order that (4.17) is tangential to the surface 7(X)=e 
we use the following conditions; 

7(X) = IxB^i - “ ^i^» j l,2,..n (4.18) 
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where is a nonzero constant. Equations (4. 18) imply that 
the slope of 'V'(Z) in the jth direction is a linear multiple 
of the slope of the hyperplane in the jth direction. 
Expanding (4.18), the following sot of equations are 
obtained . 


4- + ••• + 2PiA + ••* + 2P^n\ 


\ <^11 


2 ^ 21^1 2 ^ 22^2 + • • • + 2 ^ 21^1 + • • • + '^'^2ri-n 


4 - PqC^^^fiCcr^) + p2^22^2^^2^ + ... + c^) 


T. 

1 


0 


2i 


2^nl^l 2 Pj^ 2^2 + ••• + 2P^X^ + ... + 2P^jiX^ 


+ P-| ^ ni ^i (o'-] ) + + ... + P„G^„f„(cr^) 


'2 n 2 -^ 2'"2 


T. C . 

1 111 


m nm m' 

(4.19) 


Equations (4.19) can be compactly written in the form 

2P Z + C P f(o:) = C T (4.20) 


where p = diag (p^^), i = lj,2,...m o.nd x. is an m-vector 

m 

given by [0, 0, ... j_ • • O] • (4.19) or (4.20) together 
with (4.17) constitute a set of (n+l) nonlinear equations 
in (n+l) unknowns Z and and call for a numerical 
technique for a solution. To these equations the (m-l) 
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equations obtained from (4.4b) by deleting the ith 
equations i.e.. cr .• = C X are appended. These (n+m) 

_ I 

equations can now be cast in a vector matrix equation 
of the type 


2P -G 
T 

0^ 0 



-Cpf(a) 

a 


(4.21) 


The unknowns in (4.21) are X, and all components of o 
except which is equal to k^. let it be aosuned that 

A 

the solution of (4. 2l)yl®ld-.a a state vector X and a 
corresponding m-vector c. Substitution of these into 
equation ( 4,1.0) yields the minimum value of V(Z) on the 
hyperplane given by (4.17). An expression for this minimum 
of V(Z) explicitly in a can be obtained in the following 


way 2 

Transposing equation (4.20) and post-multiplying by 
X and reareranging, we obtain 

vn A mmA ,A. 

2X^ P X - T 0^ X = - [p f(2)] X . 


Using (4.4b) this ec^uation simplifies to 
2X PX-T a; = -[P f(2.)] £ 


Therefore 

X P X 


2 


[P f(2)]’^ a] 


Hence V(X) = rLf" 2 “ 


A 

m ^ i 

[p £(2)]^ 2 ] +■ I Pq r ^i^®i)'^®i 

i=l o 


( 4 . 22 ) 



107 


To eliminate _T consider equation, (4.20) wiiicli gives with. 

A ^ 

X = X and n = o, 

X = |-^[ax - c p f(S)] 

T 

Premultiplying "by 0 gives 

m A m "1 

0 X = P [Ct -op f(a)] 

or n = P”^[C T- C p f(o)] 

= P“^ CT - ^C^ p-^ C[p f(a)] 

. T = 2[o'^ P“^ 0]“^ o + P |.(a) (4.23) 


Substituting for 2 (4.22) V(X) is given by 

/V m 1 -1 ^ 

V(x) = /[C^ p-^ 0 ] a + I P^ J f^(Sj_)ds^ 

i=l 0 

(4.24) 


Thus V(X) is explicit in a and is useful in computing 
the values of 7(X) quickly from a knowledge of the 
sector violation information. Since there are 2m 
hyperplanes of the type (4.17) as shown in (4.16) the 

A 

minimum of all the 2m values of 'V'(X) obtained using the 
above procedure yields e which defines the region of 
attraction for the system (4.4) according to (4.11). 


4.2,3 Special cases: 

(i) The single nonlinearity case: 

The case of a system with a single 
nonlinearity (m=l) can now be considered as a particular 
(4.4). Choosing 


case of 
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A 

0 


Min(lc^, 5,^) 


equation (4.21) yields 





A 


2P 

-G 


X 


c* 

0 


T 


S 



! i 1 

ft. 


-0 P-, f(a) 


A 

a 


(4.25) 


where x and are nonzero scalars. These equations 
constitute a set of (n+l) equations in Z and x and can 


be explicitly solved to obtain ^ 

(G^P ^G) o 


(4.26) 


min' 


which is exactly identical to the results obtained by 
¥alker and Mc’Glamroch [68] and Weissenberger [67]. 


(ii) Integrals of nonlinearities neglected (p=0): 
Since = 0 (i = l,2,...m), V(X) becomes 


V(X) = P Z 

Equations (4.2l) then becomes 


2P 

-G 


— — . 


'0 

T 

G^ 

0 


X 

1 

1 

a 


(4.27) 


(4.28) 


(4.28) constitutes a set of (n-ka) linear algebraic 
equations in X, and (m-l) components of 0 that are 

A 

to be determined. Explicit expression for ■V'(X) in terms 
of the specified value for is possible in' the 


following way: 
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Por all = 0, i = l,2,...in, eqimtions (4.19) 


yield 


2P X - 0^ = 0 


(4.29) 


With. (4.17), (4.29) can be put in the form 


2P -C. 


11 


0 " 

—1 





0? 0 


I 

T . 

— 

k. 

—1 


1. ; 


1 


(4.30) 


Solution of (4.30) gives 
det[P] 

\=% 


det[H] 


^i 


(4.31) 


where H = 


2P -C^ 


h-i 


0 


and det [ • ] denotes the determinant 


of the corresponding matrix. 

If the solution vector Z of (4.30) is Z then 
simple manipula-tions in (4.29) using (4.31) yield 


2^”^ det [P] k? 


7(X) = 


det [H] 

The minimum, of all such minima oni& yields £ that defines 
the region of stability. A similar expression was obtained 
by V/eissenberger [67] for the stability domains in the case 
of a system with a single nonlinearity. 

Case (iii) Quadratic part neglected (i.e. P=0): 

In some situations it is possible to neglect the 

m 

contribution of the quadratic X P Z compared to the 
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integrals of the nonlinearity when the state variables 
assTnaed values near the point where the sector viola.tion 
occurs. In such cases it is possible to avoid the 
solution of the (n-fm) nonlinear equations of (4.2l) as 
shown below. Making P = 0, the V-function in (4.10) takes 
the form 


V(Z) = f- p. 

i=l ^ o 
and equation (4.21) reduces to 



(4.33) 


0 

-0 


"s' 


■-0 p f (a) 


0 


T 

— 

a 


from which we obtain 
G p f(a) = Or 
or C( p f (ff) -£ ) =0 
The trivial solution of (4.35) yields 
P 1 ( 2 .) =1 = [0, 0 ., .. 0] 


(4.34) 


(4.35) 


(4.36) 


Equation (4.35) implies 

f^(o.) = 0 f j = l,2,,,.m5 j T^i 

J 


(4.37) 


or 


cT, — 0 f j — J ji^i 

J 


This implies that all a.'s, (j = l,2,...ia, 2 ^±) are zero* 

V 


Also since the ith component of a is a priori specified 

in (4.17), the ith equation of (4.35) shields x j_ according 
to 

a = h 


(4.38) 



Ill 


Thus the solution vector of 0 is given by 

£ = [ 0, 0, ... k^, 0, .. 0]^ (4-39) 

and may be directly substituted in ( 4 . 24 ) to obtain a 
minimum of V(Z) on the hyparplane given by (4.17). 

Equation (4.39) thus defines precisely a point X on 
Li3> at which one of the 0 ^*s reaches its limit according 
to ( 4 . 9 ) while all others assume zero va,lues. ¥e shall 
define such points on IB as central points. 

The minimum of all the V-values computed at all 
the central points onlB thus yields e that defines the 
region of attraction according to (4.11). 

4.3 STABILITY RBG-IOFS FOR OTLTIMACHINE POWER SYSTEMS: 

4 . 3.1 Problem formula,! ion; 

In this section a method for computing the 
transient stability regions for m^ilt imac hin e power 
systems based on the results of the previous section 
is proposed. This method utilizes the sector violation 
properties of the nonlinearities of the power system 
state models. In Section 2.6 of Chapter II sta.te models 
of the fold 

Z = A X - B f ( 0 ) (4.40a) 

£ = 0^ Z (4.40b) 

have been coijstructed in the minimal state space of 
(2n-l) for the non-uniform damping and (2n-2) for the 
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uniform damping cases. Thus A, B and constitute 
a minimal rea,lization of the transfer function matrices 
W^jCs) and Wjj(s) given in equation (2.60) and (2.63) of 
Chapter II. Por such systems Lyapunov functions of the 
type 

7(Z) = P I + f ^ f^(s^)ds^ (4.41) 

have been constructed in reference [10], where P is a 
positive definite symmetric ma.trix and a constant. 

These Lyapxmov functions are used in the analysis to 
follow. The procedure of constructing such Lyapunov 
functions is based on a theorem due to Moore and 
Anderson [ 81 ]. This theorem is stated and the results 
of its application in the derivation of the Y-function 
(4.41)are shorn in Appendix 0. 

The nonlinearities i = l>2,..,m occurring 

in power system models (4.40) do not satisfy the Popov’s 
sector conditions (4.8) over the entire state space. In 
fact from (2.38) is known to violate this sector at 

= .-2(6° - a°) 

(4.42) 

and 4^ = u+2(6° - 6°) 

as shown in Figure 4,2. Therefore the Lyapunov function 
V(X) satisfies the necessary sign definite properties for 
only those cr^^, i = l,2,,..m that lie in the closed interval 
e [ii-2(6°-6°), -it-2(6°-6°)] (4.43) 

Hence a region of asymptotic stability is defined for all 
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lying in the open interval 

e {% - 2(6°-6°), -7t-2(6°-6° ) ) . (4.44) 

In this region the tine derivative V(Z) is given hy [10] 

V(X) = -Z^l X - 2 a f^(a)o; (4.45) 

for all = CO (infinite Popov sector). Ll^ is obtained 
fron the lyapnnov matrix equation 

P + P A = - 1 - (4.46) 

An estimate of this stability domain will now be obtained 
utilizing the results of the previous section. 

4.3.2 Determination of stability regions: 

The procedure of determining the stability 

region of the power system is slightly involved because 

of the fact that not all the a^’s are linearly independent. 

In fact, examination of equation (4,40b) reveals that only 

(n-l) components of the m-vector a are linearly independent. 

The remaining (m-n+l) components coni be obtained as a linear 

combination of these linearly independent, variables [24], 

Hence the boundary of the polyhedron IB defined in 

Section 4.2 is uniquely defined by these linearly 

independent variables. As an example consider a 4-machine 

system. Equation (4.40) gives 

cr^ 0 0 0 1 0 0 

02 0 0 0 0 1 0 

^ 0 0 0 0 0 1 

0 0 0 -1 1 0 
Oi- 0 0 0 -1 0 1 

Ug 0 0 0 0 -1 1 
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in whioh. the coefficient natris has only 3 linearly 

rp 

independent rows. Thus the rank of 0 is only 3. This 
inplies that only3crj^'s are linearly independent. For 
instance, ^2 '^3 assuned to be linearly 

independent. Then the other varis-bles are given by 



a. 


1 


a 


5 



and cTg 



(4.48) 


In view of this property, 
partitioned in the form 


r 1 

■ T 1 

■2i 


h 


= 




ci ■ 

-D 


D 1 

U _ 


L J 


equation (4.40b) can be 


(4.49) 


where Cj and are the linearly independent and dependent 
subvectors respectively of e* The matrices and are 
the properly partitioned submatrices of 0, The choice of 
the linearly independent set of however is not unique. 

Any (n-l) of the n variables of a whose linear combination 
results in the components of can, infact, be chosen. 
Again, as an example, consider the 4“macliine case. From 
( 4 . 47 ) one can choose a^, and as the linearly 
independent set. Then the other three variables 
and Cg are derived to be 
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(J ry (J-, CTr- 

515 

- aj_ ( 4 . 50 ) 

and cTg = 0^ - a2 + 0^ 


The numler of sucii sets of 0 j's has an upper bound of 

Each of these ffj’s thus \miquely defines a poly- 
hedron^ in the state space. Thus an upper bound on the 
possible number of polyhedrons is also 


Fow consider one such linea.rly independent Oj 
and define (2n-2) hyperplanes in the following way; 


(0^ Z)j_ = ti- 2(6°-6°) , i=l,2, ...(n-1) 

and (C^ X)^. = -iT-2(6p-6° ) , j = l,2,...(n-l) 


( 4 . 51 ) 


These equations, exactly identical to those defined in 

(4.16), define the boundary of the (n-l) dimensional 

closed surface!!^ , { Z; 'n;-2 ( 6° -6° )< (o-f) ± --n— 2( ) }., 

P 4. i P H 

This convex poljrhedronlB has (2n-2) central points* Let 
the coordina.tes of that central point at which minimiim of 

A 

V(Z) occurs be X. Then following the procedure discussed 
in case (iii) of section (4.2.3), under the assumption 
that the contribution due to quadratic part is negligible, 

A 

V(Z) is given by 


1 m 

C!j) Gj+ J fj_(Sj_)ds^ 

A [p A 

where 0j = Oj X . 


( 4 . 52 ) 
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The Eiiniinum of all such 7-values ohtainod on all 
polyhedrons yields the value of £ that defines the region 
of stability according to equation (4-. 11). Thus the 
calculation of V(Z) at a maximum of (2n-2) central 

point (120 points in a 4-machine case) is reqij.ired to 
obtain e. 

For practical systems such enumeration and 
computation is difficult when the number of machines 
becomes large . ¥e shall now discuss a procedure that 
reduces this computational complexity. One of the 
machines say ’r' is chosen at a time and the linearly 
independent variables ^ “ 1,2, . . . (n-1) , are 

selected in such a wa-y that the corresponding rotor 
angular difference (b°“6°) involve the angles of this 

ir 4. 

machine 'r', i.e., 6° or 6° = 6°. This is equivalent 

to taking machine 'r' as reference. ¥e are also assured 
that such a choice of(oj)jj^ results in a linearly depen- 
dent set Ojj. Since there are n machines in all, there 
will be 'n* polyhedrons each having its boundary JE- 
containing (2n-2) central points. Thus there will be 
a total of n(2n-2) central points at which 7(X) has 
to be computed^ Moreover, at certain of these points 
the components of obtained from the corresponding Oj 
may violate their respective sector conditions of the 
form (4.9). Such points are discarded since 7(X) may 
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not satisfy the properties of a lyapunov fnnction. The 
ninim-uia of the V-values computed at all the valid central 
points yields an estimate of e. 

The method discussed above is now summed up in 
the form of an algorithm for practical implementation 
as follows! 
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This gives a total of (2n-2) contral points. Now 
each of these central points represents one set 
of Oj. 

Step 6: Conpute for each of the central points in 
step 5 and check for violations of sector 
c conditions in the elenents of If there are 
no violations conpute V(X) , Otherwise discard the 
central point. 

Step Obtain the nininun of all V(Z) computed in 
step 6. 

Step 8 j Ghoose the next machine by incrementing in r 
(i.e. r = r+l) and go to step 4. If r=n+l 
go to step Si 

Step 9: Choose the minimum of all the minima obtained in 
step 7., This defines e. 

The algorithm was implemented on several systems 
and the results of study on a 5-machine system are 
presented in section 'Results of Study’, 

4.3.3 Comparison of the method with that of 
reference [31]* 

The technique of estimating the stability region 
detailed above will now be compared with that of 
reference [31]. T’or this to be possible, it is necessary 
to study the physical significance associated with a 
central point. As mentioned earlier, at this point on~l^ 
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all the linearly independent variables (a^) (i=l, 2. . . (n-1) ) 

i 

of excepting one assume zero value. The nonzero 

component ass-ones a value of ■k--2(5°-6°) or -ti;-2(6°-6°) . 

P q P q'^ 

For example consider the 4-machine case. The linearly 
independent set with machine 1 as reference is 

(2i)^ = 

(ap)^ = ( 6 p- 63 )-( 6 °- 6 °) (4.53) 

and (Op)^ = (6p-6^)-(6p-6°) . 

With this set one possible central point is 
(6p-62) - (bp-d^) = 'n;-2(6p-62) 

(bp-b^) - (bp-b°) = 0 ■ (4.54) 

and (bp-b^) - (bp-b°) = 0 

Consider the system on no load so that bp = 0, i=l,2,3,4. 

Then with' (4.54) one obtains bp=0 (since it is the reference)^ 
&2 - b^ = 0 and <5^'= 0. This situation iraplies that 
the value of V(X) computed at this point corresponds to the 
potential energy required for talcing the machine 2 out of 
step from the rest of the system. Thus the system while 
losing synchronism is divided into two groups one of which 
contains machine 2 and other contains the rest [l7»3l]. 

This analysis applies to every other point. It is 
precisely on this basis that the authors of reference [3l] 
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approxinated the xmstable equilibriaiTi points. It nay 
therefore be concluded that the method proposed in 
Section 4,^.2 provides a necessary theoretical rationale 
for the development in reference [3l] with the following 
additional features; 

1. The search among the n(2n-2) central points of 'n' 
polyhedrons takes into account the machine reference 
automatically. 

2. An estimate of e is obtained without violations of 
the sector conditions of any of the n nonlinearities . 

The method, however, is found to yield slightly 
conservative reailts compared to those of reference [31]. 

4. 4.3.4 An alternate procedure of estimating the 

region of stability; 

An alternate procedure utilizing the sector 
violation properties of the power, system nonlinearities 
that yields better stability domains is suggested here. 
Considering the fact that only the first (n-l) Variables 
of a with machine 1 as reference are linearly independent, 
it was conjectured in reference [24], that the minimum 
of occurs at one of the vertices of an (n-l) 

dimensional polyhedron in the cr-space defined by this 
set of linearly independent variables accoiding to 
equation (4.51). But the fact that any (n-l) components 
of £ can be chosen to be linearly independent as described 
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earlier was overlooked* Using ihis proijerty it is now 
possible to estinate the region of stability as follows s 

Consider a polyhedron defined x^ith nachine 'r' as 
reference, i.e. either 6° or 6° = 6° in the set (4.51)* 

This polyhedron has 2^“^ vertices. It can be shown that 
each of these points correspond to the nachine 'r' going 
out of step with respect to the rest of the machines. For 
example, consider equations (4.55) of the 4-nachine case 
where machine 1 is the reference. The polyhedrorUB formed 
by the 5 linearly indepdndent variables (o;j)^ » i=l,2,3 ha,s 
8 vertices. Also assume an unloaded case such that 6? = 0, 
i = 1,2,5»4* Then the coordinates of the vertices of are 
(-n:. It, It), (-jt, -Tt, Tt) and so on. Bach of these vertices 
therefore, corresponds to the machines 2,3 and 4 going out 
of step with respect to nachine 1. 'Equivalently nachine 1 
looses synchronism. Moreover, the linearly dependent 
variables can be verified to be taking value of 2it or 0. 

This implies that the dependent variables are pinned down 
at the origin of their respective cr~f(a) planes since the 
solutions are periodic in 2Tt, 

In a loaded system, however, the minimum of the 
values of V(X) computed at the 8 vertices oflB yields the 
minimum energy required to dissocia,te machine V from 
the system. Similar analysis holds for the other machines also . 
Since there are n machines in all, the minimum of all the 
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7-values obtained on all tlie polyhedrons yields an 
ostinate of s tiia,t defines the region of stability. 

In view of the above analysis, the following 
algorithm is proposed. 

AI&ORITHM 2; 

Step 1: Obtain the sta.ble equilibrium point of the 

post-fo,ult system 2,3? •••n. 

Step 2: Calculate all possible angular differences 

(6^-6°), p,q = l,2,...n (p< q) . The total number 
P Q, 

of such differences is n (= n{n-l)/2). 

Step 3' Choose machine ’r'. Set with r = 1, It is assumed 
here that the machine 1 goes out of step. 

Step Among the differences of stop 2, pick all the 

(n-l) differences such that p or q = r. 

Step Compute the values of k^ and for each of 
the (n-l) differences in step 4. 

Step 6: Choose 

= if Si^< k^ , i=l,2, . ..(n-l) 
and form the vector Op. 

Step Jt Set all the (n-n+1) linearly dependent variables 
at -2(6°-6°), and form the vector 
step 8: Compute 7(X) . Set r=r+l and go to step 4 if all 
the machines are not exhausted. If r = n+1 go 
to step 9. 



124 


Step 9: Choose the nininm of all V-values obtained in 
step 8. This gives e. 

This nethod was inpleaented on a nmiber of systems 
and has been foimd to give highly accurate results. Further, 
the nethod avoids the problem of machine reference as in 
method I and the value of 'V(x) for a particular ’sa chine 
loosing synchronism is uniquely determined. Results of 
study on a 5 -machine system along with a. comparison with 
the other methods are given in the next section. 

4.4 RESULTS OF STUDY: 

The system chosen for study is a 5-machine system 
the data of which is given in Appendix D, The magnitudes 
of the voltages behind the transient reactances, mechanical 
power inputs and inertia constants are given in Table 4.1. 

A fault at bus 6 is cleared by isolating line 15. The 


Table 4.1: System Data 


Machine 

Voltage behind 
transient react. 

(Ej_) 

Mech. power 
input 

'•"mi^ 

Inertia 

constant 

(%) 

1 

1.05754 

0.900 

2.50 

2 

1.08650 

1.750 

10.50 

5 

1.08521 

5.200 

12,72 

4 

1.05586 

^ 5.500 

27.00 

5 

0.95777 

0.445 

oo 
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driving and transfer admittances between the internal 
nodes of the machines are calculated and shovn in 
Tables 4.2 and 4.5 for the fa'olted and the post-fault 
systems respectively. Transfer conductances are neglected 
for the post-fault state but the short circuit conductances 
are read^justed to satisfy eQ.uation (2.14) of Chapter II. 

All data are in p.u. 


Table 4.2(a) : Gr-matrix for the faulted system. 


M/G Ho. 

1 

2 

5 

4 

5 

1 

1.17664 

- 0.00560 

0.0000 

-0.07078 

- 0.50952 

2 

-0.00560 

0.25677 

0.0000 

0.08580 

- 0.00924 

5 

1 

0.00000 

0.00000 

0.0000 

0.00000 

0.00000 

4 

-0.07078 

0.08580 

0.0000 

5.42242 

-0.27819 

5 

-0.50952 

- 0.00924 

^ 0.0000 

-0.27819 

1 2.56496 


Table 4.2(b): B-matrix for the faulted system. 


M/C Ho. 

1 

2 

5 

4 

5 

1 

-4.52778 

0.28854 

0.0000 

-0.94978 

1.86827 

2 

0.28854 

- 5.42500 

0.0000 

0.80560 

0.50746 

5 

0.00000 

0.00000 

- 15.9470 

0.00000 

0.00000 

4 

.0.94978 

0.80560 

1 

0.0000 

-10.88625 

5.88867 

5 

1.86827 

0.50746 

0.0000 

5 .88867 

-8.80014 
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Table 4.3(a)- G-matrix for the post-fault system 


M/0 Mo. 

1 

2 

3 

4 

I 

5 

1 

1.17877 

0.01445 

0.04455 

-0.03210 

-0.50658 

2 

0.01445 

0.29015 

0.15703 

0.16609 

0.02151 

3. 

0.04455 

0.15703; 

1.33232 

0.35635 

0.07116 

4 

-0.03210 

0.16609 

0.35635 

3.62085 

-0.21297 

i 

5 1 

-0.50658 i 

0.02151 

0.07116 

-0.21297 

2.36799 

I 


Table 4.3(b): B-matrix for the post-fault system 


M/0 Mo. 

1 

2 

3 

4 

5 

1 

-4.42707 

0.38725 

0.71489 

1.22917 

2.04542 

2 

0.38725 

-3.33079 

0.69596 

1.07222 

0.68124 

3 

0.71489 

0.69596 

-5.16262 

1.97359 

; 1.25780 

4 

1.22917 

1.07222 

1.97359 

-10.12185 

4.38045 

5 

2.04542 

0.68124 

I 1.25780 

4.38045 

-8.56853 

j 
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The post-fault eciullihritua state (with the transfer 
conductances neglected) is given by 


c ^ c O 

«l-*2 = 

-0.34147 

C O c 0 

-0.19861 

aj^-64 = 

0.04049 

xO A® 
6^-6^ - 

0.08604 . 


Damping is assumed to be uniform and X is chosen to be 
unity. Then the swing equations describing the system in 
the post-fault state with machine 1 as reference are given 
by 


di 

dt^ 




5 

- I sin(6.-6.)}], j=2, 3,4,5 (4.55) 

i=l J J 


The parameters in (4.55) are obtained from Tables 4.1 
and 4.3. When these equations are cast in the Lure’ form, 
the number of nonlinearities f;j^(o’j_) is 10. They are 

fl(cri) = {sin(cr^ + ( 6 ^- 6 ^)) - sin( 6 ^- 62 )] 

^2^ ^2^ = [sin(o2 + (6^-6^) ~ sin(6j-6^)] 

f^(o-^) = [sin(cr^ 4 ( 6 ^- 6 ^)) - sin( 6 ^- 6 ^)] 
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= [slnia^ + (6j-6°)) - aiii(6j-6°)] 

= tsin(a^ + (62-6^)) - sirL(62-6°)] 

^6 ^^6^ = [sin(ffg + (62-6°)) - 8111(62-6°)] 
f^Cc?^) = [sin(a^ + (6°-6°)) - sin(6°-6°)] 

= [sin(crg + (6°-6p) - sin(6°-6p] 

^9(^9) = [sin(cyg + (53-6°)) _ sin(60-6°)] 

^ [sin(<j^Q+(6°-6°))- sin(6j-6^)] . 

These nonlinearities violate their corresponding sectors 
at and — given in Table 4.4. 


Table 4.4- Sector violation information 


SI .Ho . 

^i 

Value of 

(6j-a°) 

^i 

ir-2(67-6p 

-2,j_ - 

-TT - 2(6° - 6°) 

J- J 

1 

^1 

-0.34147 

3.48306 

-2.80012 

2 

02 

-0.19861 

3.34020 

-2.94298 

3 

0^ 

0.04049 

3.10110 

-3.18238 

4 

<^4 

0.08604 

3,05555 

-3.22763 

5 

^5 

0.14286 

2.99873 

-3.28445 

6 

^6 

0,38196 

2.75963 

1 -3.52555 

7 

Oy 

0.42751 

2.71408 

-3.56910 

8 

<"8 

0.23910 

' 2.90249 

-3.38069 

9 

^9 

0,28465 

2.85694 

-3.42624 

10 i 

°^id 

0.04555 

3.09604 

-3.18714 
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Here four of the ten nonlinearitieo are linearly 
independent. Algorithm 1 will now be implemented. 

Choosing the machine 1 (r = 1) first, the linearly 
independent variables selected are o-2_» ^ 2 * ^3 ^4-* 

With these as components of Oj, one central point is 
obtained with the following values ^ 

Oj = -2.80012; 02 = 0 ; = 0 and cr^ = 0. 

This point corresponds to machine 2 going out of step 
and machine 1 is the reference. With this set const i~ 
tuting Oj, the components of are : 

= 2.80012 ; Cg = 2.80012 j = 2.80012 

cTg = 0.0 ; = 0.0 and = 0.0 . 

With these, V(2) is calculated to be = 3.42169. 

It may, however, be noted that Og and Oj violate 
their sector conditions. According to the algorithm, 
this point is inadmissible. Table 4.5 shows such an 
analysis at the other centra.1 points. Only 4 points 
have been shown because the rest four points on this 
polyhedron yield larger values. Similar procedure is 
adopted for other machines ’r’, r = 2,3»4 and 5. A 
summary of this analysis is shown in Table 4*6 where the 
values of V(X) for the 5 machines going out of step and 
with every machine selected as reference are given. Also 
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the sector violation conditions a,ro indicat od. The 
miniiiinin of all the V-values at all the admissible 
central points is 3.18156. Hones the stability boundary 
is defined by 


V(X) = 3.18156 . 

Table 4.5* Choice of M/c 1j An analysis of central points. 



Now the procedure indicated in algorithm 2 will be 
adopted. Accordingly we first choose machine X (r=l). 


Hence the components 

are assigned values 

of Cj are cr 2 _» ^ 2 * ^3 

= -2.80012 ; 

= 3.1011 

^2 = -2.94980 ; 

<7^ = 3.05555 . 

Assign values of -2(6 

?-6^) for all the other linearly 

dependent variables. 

Thus 

= -0.28572 J 

ag = -0.47820 

ag = -0.76392 ; 

Og = -0.56930 

Orj = -0.85602 ; 

010= -0.09110 
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Table 4.6: Effect of nacliine references. 


SI. 1 
ETo. i 

1 

Hef . 
M/c 

Variables 

of Oj 

EMachines loosing synchronism 

1 2 3 4 5 

1 



-2.80012 0.0 0.0 0.0 

2 

1 

^2 

0.0 -2.94298 0.0 0.0 

3 


0.0 0.0 3.10110 0.0 

4 


^4 

0.0 0.0 0.0 3.05555 


Value 

of V(Z) 

3.42169 7.50162 15.71059 13.87118 



(VL) (VL) (VL) (VL) 

5 



-2.80012 - 0,0 0,0 0.0 

6 

2 

^5 

0.0 - 2.99873 0.0 0.0 

7 


^6 

0.0 - 0.0 2.75963 0,0 

8 


arj 

0.0 - 0..0 0.0 2.71408 


Value 

of V(Z) 

8.75572 - 12.8044 15,58135 13.73607 



(vn) (vn) (VET) (vn) 

9 


''2 

-2,94298 0.0 - 0.0 0.0 

10 


"^5 

0.0 2.99873 - 0.0 0.0 

11 



0.0 0.0 - 2.90249 0.0 

12 


'Oq 

0.0 0.0 - . 0.0 2.85694 


Value of V(Z) 

8.90831 3.18156 - 15.7627 13.89756! 

(VL) (VET) (VL) (VL) 

13 


^3 

3.1011 C.O 0.'^ - 0.0 

14 

4 


0.0 2.75963 0.0 - 0.0 

15 


O.Q 0.0 2.90249 - 0.0 

16 


^10 

0.0 C.O 0.0 - 3.09604 


Value of V(Z) 

9.46432 3.45545 7.52921 - 13.81387 

(VL) (VL) (VL) (VL) 

17 


^4 

3.05555 0.0 0.0 0.0 

18 

5 

4 

0.0 2.71408 0.0 0.0 

19 


0*9 

0.0 0.0 2.85694-0.0 

20 


°’io 

,0.0 0.0 0.0 3.09604 


Value of V(X) 

9.44373 3.48729 7.55031 2 0 B99p - 

(VET) (VET) (VET) (VET) 


note: VL : indicates sector violations-. 

VET : indicates no sector violations. 
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With these variables V(Z) = 4.03076. 

The procedTire is followed in a similar manner for the 
other machines i.e., r = 2, 3, 4, and 5. The Y— values obtained 
at the corresponding corner point are given in Table 4.7. 
The minimum of the V-values o.t these points yields e. 


Table 4.7 J Comer point evaluation. 


SI. M/c going linearly independent variables Value of 

Wo, out of step to be chosen V(X) 


a-. = -2.80012 ; cr, = 3.10110 I 
1 1 ^ 9.68689 

0 = -2.94298 ; = 3.05555 


cr, = -2.80012 ; a. = 2.75963 I 
2 2 -^ ^ I 4.03076 

CT5 = 2.99873 ; Orj = 2.71408 


a = -2.94298 j 0« = 2.90249 
5 3 ^ ® 7.66021 

= 2.99873 ; = 2.85694 


cr~ = 3.10110 ; Oo = 2.90249 

4 . ^ 15.60612 

erg = 2.75963 ; cy^Q= 3.09604 


a. = 3.05555 ; o-o = 2.90249 

5 5 ’’ ® 14.84772 

Orj = 2,71408 ; 03^3= 3,09604 


A comparative study of the various methods 
including that of reference [3l] will now be made. 
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Table 4.8 gives the values of the Y-fxmction and the 
corresponding estimate of the critical clearing times 
obtained. The actual critica.1 clearing time is also 
obtained by integrating the differential equations of the 
■fcype ('’-.55) in the faulted state. The parameters in these 
equations assume . values given in Tables 4.1. and 4.2. The 
critical clearing time computed, is also shown in Table 4.8. 


Table 4.8: A comparative study. 



Values of V(X) computed using method of 


ref , [ 5l]* 

Algorithm 1 



Algorithm 2 

Actual 

integration 

Value of 

v(x) 

5.48729 

' 

5.18156 

4.05076 

- 

Critical 
clearing 
time sec, . 

0.1780 

0.175 

0.190 

0.205 


* Choosing the machine with the largest inertia as reference. 


It is observed that the critical clearing time 
computed using the results of algorithm 2, yields clearing 
times close to the actual, than those offered by the 
technique of reference [51] and the algorithm 1 s^uggested 
in this chapter, 

4,5 COICIUSIONS: 

In this chapter regions of attraction for multi- 
nonlinear systems using the sector violation properties 
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of tile nonlinearities have been derived. These results 
constitute a generalization of those obtained in the 
literature for systems with a single nonlinearity. This 
generalization involves the solution of. a set of nonlinear 
equations in a general case. It is also sho^^m that explicit 
expression is possible in special cases. These results 
are extended to the multimachine power system stability 
problem and it has been shown that they provide a theore- 
tical basis for a method suggested by Prabhakara and 
El-Abiad [3l]. However, this method is suitably modified 
by considering the, sector violations of the dependent 
nonlinearities and machine reference. The procedure is 

summarized in algorithm 1. Algorithm 2 is new. Both of 
these have been implemented on a number of practical 
systems [ 79]. Results of study on a 5-maGhine system 
are given and the methods are compared. 



CHAPTER Y 


COHERBHCT IDEHTIPIGATIOH AM) DYHAMIO EQUIVALENTS 
USING ENERGY FUl^CTIONS 


5.1 INTRODUCTION! 

This chapter deals with the development of reduced 
order models of power systems for use in the transient 
stability of large scale power systems. Due to the 
increased interconnections in power systems, it usually 
happens that the behaviour of one section of the system 
may not affect significantly the performance of some other 
section of the system under study. The stability analysis 
of this latter section called the 'study system' can then 
be carried out by a simplified representation of the former 
section called the 'external area' both in size and in 
detail. This proced-ore reduces the burden on the computer 
memory and saves considerable computation time. The reduced 
order model of the external system thus obtained is called 
the 'Dynamic equivalent'. 

The development of dynamic equivalents for large 
power systems for use in stability analysis, however, is 
not as easy as for load flow studies [80], Some of the 
difficulties associated with equivalencing are as follows 
[ 55 ]- 

(i) The dynamical nature of the synchronous machines, 
governors etc, cannot be easily represented on the 
equivalent. 
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(ii) A proper identification of the portion of the 
system to be equivalenbed is required, 
and (iii) Inertia allocation among the equivalenced machines 
for multimachine representation poses formidable 
problems. 

Currently, three different approaches exist for 
obtaining the dynamic equivalent for a large power system. 
These are 

(i) Use of netvrork; equivalents and inertia alloca- 
tion [5'4-, 55 ], 

(ii) System reduction employing the property of 
coherency among machines[55»56 ,57,59 ], 

and (iii) Modal techniques [61-64]. 

Of these, coherency based dynamic equivalents are by far 
the simplest. This has a potential application as a 
production tool in planning and represents a key step in 
the on line monitoring techniques [57]. The method 
consists in identifying the machines that swiug together 
during a transient disturbance by processing the swing 
curves. Pairs of generators outside the study area that 
tend to swing together are compared and generators are 
classified as coherent if the maximum excursion in their 
angular difference does not exceed a pre-assigned limit[55]. 
This concept of coherency has been extended to generator 
as well as non— generator buses of the power system [56] » 
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Two buses are defined as coherent if for all time (i) the 
ratio, .of voltage magnitudes of the buses is constant and 
(ii) the difference in voltage angles of the two buses 
is constant. A method for coherency identification without 
actually conducting a stabilitjr study has been suggested 
by Lee and Schweppe [59]. They utilize the concept of 
’a distance measure' to ascertain the importance of a 
generator to the stability of the machines in the study 
area and propose a Set of 'features' a term used in 
pattern recognition, for identifying coherent machines. 

Once coherency among the machines is established, there 
remains the problem of deriving the dynamic equivalent. 
Chang and Adibi [53] have used the well known node elimi- 
nation formula. This technique alters every element of 
the admittance matrix with all the internal nodes of 
generators retained. A simple computational procedure 
based on the power invariance criterion has been proposed 
by LeMello etal [56,57]. Here only those elements of the 
admittance matrix that correspond to the coherent group 
are altered. 

In this chapter energy function is proposed as 
a criterion for coherency recognition as well as for 
deriving the dynamic equivalent. These energy functions . 
have been used as Lyapunov functions in the transient 
stability investigations of power systems. The reduced 
order dynamic equivalent derived by using this criterion 
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is iisefuL in the stability analysis of the imreduced system 
via the Lyapunov methods. There are also advantages of such 
an analysis in the evaluation of the unstable equilibrium 
points of the post-fault state of the system as pointed out 
in a- recent paper by luders [ 17 ]. ■ 

A useful step in obtaining the dynamic equivalent 
is to divide the system into the following two areas: 

(i) the faulted area or usually referred to as the 
study area in which the system behaviour is of 
direct interest. This is usually specified 
a priori, 

and (ii) the external area or the remainder of the system. 

The dynamic equivalent is derived for the external area for 
disturbances in the study area. This determination is 
carried out in two steps as shown in figure 5.1. These 
are (i) coherency recognition and (ii) dynamic equivalencing 
and will be outlined in detail in the following sections. 

5.2 GOHEHEICY REOOGhITION: 

5.2.1 Conventional definition of coherency: 

Traditionally, the term coherency has been defined 
with respect to machines that tend to swing together. Two 
generators ‘i’ and ’j* are said to be coherent [ 53 ] if 
there exists a constant such that 

6i(t) - 6j(t) = ajj 


(5.1) 
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for all tine. Similarly a group of generators is said to 
be coherent if every pair of machines in the group is 
coherent, for example, consider a ss^stem consisting of 
6-machines whose swing curves are as shown in Figure 5.2. 

It is evident that the angular difference 62_2 t>ofween the 
machines 1 and 2 is constant subseciiient to the disturbance 
at all time, i.e., 

6i2(to) = 

NOTE: In this chapter double subscript notation is 
introduced to denote the angular differences between two 
machines . 

Similarly machines 3 Q-ad 5 maintain a constant 
angijlar difference. Thus two coherent groups with 
machines 1 and 2 in one and 3 and 5 in the other can be 
formed . 

Differentiating (5.1) with respect to ’t’ and 
roar ranging the terms one obtains 

= “("fc) (say). (5.3) 

Therefore if two or more machines are coherent, then their 
a,ngular velocities are the same. The kinetic energies 
associated with the machines i and j are 

K.E.j_(t) = i |(t) 

E.E.j(t) = i w2(t) 


(5.4) 
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Also the potential onorgy P,E.j_j(t) associated with 

this pair is given by 

a. . 


' 0 - 13 ^ 13 -13 13' 

-S 


- oos{e^ -6..)] dOj 


13 


( 5 . 5 ) 


where Y. ./Q- • = 1- • + 3B. - is the short circait transfer 

13^ 10 10 10 

s 

admittance between the machines i and 3. 6^^ is the 

pre— fault stable angular difference between the two 
machines i and 3. is given by 




( 5 . 6 ) 


Under the assumption that all transfer conductances are 

neglected (5.5) reduces to 

a. . 

P.E.ij(t) = J B^B.B.j[sin(ti^j+6^p-sin 6=^] da3_j 

7 ’> 


The total energy ^(t) associated with this pair of machines 
at any tine 'f is then given by 




•^(t) 


i M.a)| 


+ 


M.Mt 
J J 


‘’ij 




- sin . (5.8) 


Equation ( 5 . 8 ) thus gives the total energy gained by this 
pair of machines i and 3 as the system evolves d'oring a 
transient period fron its pre— fault state. If however. 
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replaced by 6?^, the post-fault steady sta.te 

angular difference between the nachines i and j we obtain 

cr. . 

- sin 6 A] . (5. -9) 

o ^ 

where o . - = 6^ ~ . . and B. . is the transfer admittance 

between machines i and o in the post -fault state. This 
energy function in (5.9) has been extensively used as a 
Lyapunov function V(X) in the stability analysis of power 
systems [5,6] and represents the energy gained by the 
system with the post-fault equilibrium as the datum. 
Either of the two equations (5.8) or (5.9) can be used 
for coherency recognition. In this chapter, however, 
equation (5.8) is used in wha,t follows. 

When the machines are coherent at all time, then 

5.Ai;) = 6®. for all t > t (5.10) 

X J X J w 

and will also be equal to 6?- when a post-fault stable 
equilibrium state exists. Therefore the integral in 
(5.8) reduces to zero and the energy '^(t) is given by 

■^(t) = i M^u|(t) +iM^.a)2(t) (5.11) 

and in view of (5.3), 

+ M.) 0)2(t) . (5.12) 

It is thus purely kinetic in nature. Hence the coherent 
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group is equivalent to a single machine whose inertia is 
the sum of the inertias of the individual generators. This 
is indeed the conventional method of obtaining the 
equivalent inertia of a single machine replacing a 
coherent group. 

The fact that the integral vanishes provides an 
alternative criterion for identifying coherent machines 
which may be stated as follows; 

5.2,2 An alternative criterion for coherency: 

Two machines i and 3 are said to be coherent if 
the potential energy component associated with this pair 
of machines is zero for all time. Mathematically 
"^i j 

° ( 5 . 13 ) 

Similarly a group of machines is said to be coherent if the 
potential energy component associated with every possible 
pair of machines in the group is zero. Thus for all possible 
angular differences 6 ^^(t) for t> t^, ( 5 . 13 ) is computed at 
every time step of integration of the differential equations 
and compared with a pre-assigned value. 

During the computation it is possible to obtain 
extremely low values for the left hand side of ( 5 . 13 ) if 
the value of B. . is extremely small. Such a situation 

1 J 

occurs when the two machines under consideration are 



electrically very far apart. To serve as a check whether 
or not the two machines are coherent , the slopes of the 
swing curves which yield the velocities are also monitored 
and compared. Instead of directly comparing the frequencies, 
the total kinetic energy is compared with the kinetic energy 
obtained by averaging the velocities, i.e. the expression 

i M,4(t) + i M.m^Ct) - (5.14) 

11 J J J 

where w (t ) = t[ (t) ] 

is computed. This is equivalent to comparing the kinetic 
energy of the equivalent machine with the sum of the 
kinetic energies of the individual machines. 

Pei’fect coherency in a practical system is very 
rare. It is therefore necessary to specify admissible 
tolerances on both the criteria (5.13) and (5.14) as 
follows : 

P.B.j_^(t)< (5.15a) 

[K.B.,,(t) - (K.E,.(t) + .^.E..(t)] < e (5.15b) 

where s-j^ and are the pre—assigned tolerances. Coherency 
identification using the criterion (5.15) is now possible 
as the stability study proceeds. At the end of every 
integration step the expression on the left side of 
equations (5.15) are computed for all possible pairs of 
There will thus be n(n-l)/2 possible 


machines . 
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computations to be carried out and those pairs that satisfy 
(5.15) are marked as coherent. ?rom the sets so obtained, 
coherent groups are formed. As an oxample, consider a 
5-ma.chine system. If the criteria are satisfied by the 
pairs (2,3), (2,4) and (3,4) at a particular time step, 
then the coherent group of machines 2,3 and 4 is formed. 

At the end of the study, all groups that appear as 
coherent in every time step of integrevtion are chosen. 

The procedure was implemented on a 15 machine, 44 bus and 
56 line system and the results are summarised in 
Section 5.4. 

5.3 DEVEIOPMENT OP THE DYMAMIC EQUIVAIBNT s 

In this section expressions for the admittance 
parameters in the dynamic equivalent obtained by replacing 
the coherent group of machines in the external area by an 
equivalent generator is obtained. Although results obtained 
for these parameters are identical to those of De Mello 
et al [56] the method of deriving the eqTxi valent is new 
and utilizes the energy criteria developed in the previous 
section for coherency identification. 

5.3.1 System descriptions 

With the usual assumptions of (i) a constant voltage 
behind the transient reactance and (ii) all loads represen- 
ted by constant impedances to ground, the self and transfer 
admittances between the generator internal nodes are first 
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derived. These are obtained bj first angmonting an 
additional set of nodes to the bus admittance matrix of 
the system. These additione.l nodes represent the internal 
nodes of the generators. Then all the nodes excepting the 
generator interno.1 nodes are eliminated. This results in 
the relationship 

I = Y I (5.16) 

where I is a vector of complex current injections at the 
generator internal nodes and B is a vector of complex 
internal voltages. A (-) over the variables represents 
complex quantities.. 

When the power system consists of groups of coherent 
generators it is possible to reduce the system model (5.16) 
by replacing es^ch of these coherent groups by an equivalent 
generator. Bor example assume that in the 4-machine ca.se 
of Bigure 5.3(a) machines 3 and 4 form a coherent group. It 
will first be observed that a boundary between the coherent 
and the non-coherent areas through the nodes 3 and 4 can be 
drawn. This divides the system into two parts, one 
containing the non-coherent machines 1 and 2 and the other 
containing the coherent groups (3>4). The generator nodes 
3 and 4 corresponding to the latter group lie on this 
boundary as shown in Bigure 5.3(b) and can be replaced by 
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a single node 'e' and the corresponding generators by a 
single generator G as shown in Figure 5.3(c). 


Assuming that only one coherent group of 'r' 
machines exists while the others swing independently 
it is possible to recast (5.16) in the form 



1 

1 

IhI 

1 



where 


% 




E 


^OC 


Y T 

m 


r 


i Sc 


E, 


-E 


(5.17) 


•¥ 


= vector of complex current injections at 
the boundary or coherent internal nodes 
and of order r. 

= vector of complex current injections at the 
rest of the internal nodes and of the order 
m( = n-r) 

= vector of complex voltages at the boundary 
nodes 

= vector of complex voltages at the non- 
coherent nodes. 


Subscripts 0 and N refer to the coherent a,nd the non- 
coherent groups respectively. In the example of Figure 
5.5(a) r corresponds to machines 3 and 4. 

If the coherent nodes are to be replaced by an 
equivalent node e, the system representation will be of 


the form 
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ielT 



(5.18) 

% 


Y 

-le 

! 

i 




wliBZ'B I ~ cixrx*6iiij in j Gcti-OH iiiG GQ_niv 0.16111; nodo 6 
0 

B = voltage a,t tlie node 'e'.. 

0 

Comparing (5.17) and (5.18) it may be noted that the terms 

in unchanged. Explicit expressions for Y and 

Y.t will now be obtained. As will be sho^m later, Y is 
— Ue 

T 

made etiual to . 

5.3.2 Computation of Y^^ and 


Expanding (5.17) one obtains 


Hi I 

a 

H 

J 


^C1,C1 '^C1,C2 **• ^Ca,Cr i^01,lTl '^Cl,lTm 

®ai 

^C2 


^C2,C1 ^02, C2 *" '^C2,Cr|^C2,Kl •** ^C2,Nm 

Eq 2 

•- 

• 

^Cr 


* * * ' „ L 

^Cr,Cl ^^Cr,C2 *** '^Cr,Cr,^r,m 

®Gr 

^E1 


%1,C1 %1,C2 **• %l,Gr "^171, N1 **•'%!, -Hm 

\l 



%2,C1 %2,C2 ••• %2,Cr‘%2,Nl •** %2,Em 

%2 

J 


%m,Cl '^IIm,C2 *** '^llm,Cr! %m,Nl %m,lj 

' 

But 1 

Hm 


( 5 . 19 ) 


Subscripted numbers in (5.19) denote the nodes in the 
respective groups U or C. The injected currents at 
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any ith maciiine in group 0 and ai any iiiii node in group N 
are respectively given "by the following eq_uations. 


^Oi = ^Ci,Gl ^G1 ^Gi,G2 ^^^ + • • • + ^Gi,Ci ®Gr ^ 


^Gi^ta %1 + ••• + ^Gi,Km 


m 


^Ci,Ga ®Gd ^Ci,Nk V 


( 5 . 20 ) 


and 


m 


E, 


y Y „ '®G1 -“Uk 


( 5 . 21 ) 


From equation (5.18) for the reduced system the current 
injection at the eth node is given by 


m 


= ^ee ^e,Nk %k 

and in the group I, at tne ilth node 

m 


(5.22) 


It.t 0 = f Tsrs, ^ 




V=1 


IfJi.Jfk %k 


(5.23) 


¥e now seek to develop expressions for Y^^ and 

Y utilizing the criterion developed in Section 5.2.2 
— Ee 

for coherency which states that the system energy in the 
reduced and the unreduced system be the same. The total 
energy "^(t) of the unreduced system at any time t with 
transfer conductances neglected is given by 
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Da +r m +r 

'^(t) = I _ i I + l_ 


i=l 


m-tr-1 m+r ij o 

i=l J=i+1 o D. J IJ J-J 


sin 6j_j]aaj^j • 


(5.24) 


( 5 . 24 ) may be written, as 


m r 

g(t) = i i 

ielT keC 

r 1 r ‘ 

+ J Ej^l!jB^.[sin(a. j+6|3)-sin6"^j 

ieC jeC 


m-1 m 


'kil 


+ 1 , , 1 , ! 

k=l fi,=k+l o 

keir £eU 


S, 


a. . 

r m ij 


P =± O 

peO q£N (5.25) 


It may be noted that the subscripts N and 0 for the 
voltages, denoting the groups to which these voltages 
belong as sho-^jn in (5.19) are omitted for convenience. 

The first two summation terms in (5.25) represent the 
kinetic energies of the machines in the groups K and 0, 

The third term is the potential energy in the coherent 
group C and is therefore zero. The foiirth is similar 
to the third but pertains to group H but is not zero. The 
last term vrhich is of interest represents the potential 
energy associated with every, pair of machines obtained by 
from C and the other from N. 


choosing one 
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In a similar manner for the equivalent, the energy 


-|*(t) is given hy 


-IMt) = 


i I M-co^t) + t M 
i=l ^ ^ ® 


’Itl 




m-l m kS, q s 

+ 1^=3, i 

kell ££11 
a 

m eq s q -, 

‘"1=1 ®eVe<lf®“'®e4-«eq>-=^"^^e<i^‘l'’eq 


— eq-^'‘‘''eq 


(5.26) 


The first and the third terms in (5.26) correspond to 
machines in N and are the same as the first and the 
fourth terms respectively in (5.25). We note that the 
coherent group can "be replaced hy an equivalent machine 
of inertia M given hy 


^e = ^ ^^i 

® i=l 


(5.27) 


and also Uv(’*")» (keC) = co(t). 


Now equating the energies ^*(t ) to^(t) we get 

O' 

qll ®eVeq["^'''°eq-«eq)-=“®Jq 

qsN ^ 

r m pq s \ s t 

peC q£N (5.28) 
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Por any arbibrary q in N, aquation (5.28) yields 
a 

eq 

f BE B [sin(a +6 ^)-sin6^ ]da 
•' eqeq^ an rh' en-^ fin 


eq eq 


)q eq 


^ PI q c. 

” pll 0^ = • <5-29) 

peG 

B is not Icnown a’ priori. It is arbitrary and is chosen 

as some sort of an average of the voltage at the buses in 

the coherent group. Different methods of selecting 

have been proposed in reference [56]. In this analysis, 

the following choice of E overcomes the difficulty in 

solving for from (5.29)* 




l6l 

e e ^ e 


pii p/ P^ 


(5.30) 


Equation (5.30) implies that E^ is the geometric mean 
of the complex voltages of the coherent nodes. With this 


choice it is trivial to show that 

a = a , P = 1.2, ...r. (5.3l) 

pq eq^ ^ 

Hence equating the integrands on either sides of (5.29) 
one obtains 


® sln( c-eq+«ea> 




(5.32) 


Since and 


one obtains 
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from (5.52) 


1 B 

sin 6^ - E E B 

q eq 

eq e q eq 

r 

r 

Y B^ 

E B^ sin 6 - 7 ; 

^=1 P 

q pq pq pii 

pec 

peC 


The first tem on the left side of (5.35) represents the 
po¥sr transfer from the eq.uivalent bus to the node' q' in 
ST and the second teim represents the power transfer at 
the initial time t^. Similarly the first and second 
terms on the right side of (5.33) yields the power 
transfer from the coherent group to node ’q' at 't' 
and 't^' respectively. Since energy transfer is to be 
same at all time, the power transfer at all times should 
also be equal. Hence the first terms on either sides can 
be equated to give 


B^E sin 6^ = I B^ B^ sin 6^ 

e q eq eq P q Pq Pq 

pea 

or Re[l f * Y* ] = \ Re[E E* Y *] 

a ® ^ P q.P 


(5.34) 


HOTB; * indicates conjugation. 


Therefore, 


q e qe 


7 1 * 1 ? 
pil q P qp 


from which 
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f = Y (E /S ) Y 

1© pil ^ pq 


(5.35) 


Thus the qth element of is obtained. This result 

agrees with that obtained by Bree etal [56]. It now 
remains to obtain ? 

ee 

To compute Y it will be assumed that the net 
energy fed into the unreduced system by the coherent group 
is exactly equal to the energy fed into the reduced system 
by the equivalent machine at all time. This implies that 
the net power injection 'P ' at the equivalent node 'e* 
in the djrnamic equivalent equals the sum of the net 
injections P^, i = l,2,...r at the internal nodes of the 
coherent machines in the unreduced system. 

The procedure therefore becomes identical to that 
by Bree etal [56]. 


The power P^ injected at any node 'i* in 0 is 


given by 


^i “ ®Ci ^Ci ' 


(5.36) 


The total power P injected into the system by 0 is 

o 


.L 1*01 ®oi - 

1=1 


(5.37) 


and the power injection P at the equivalent node is 

© 


p = r * 

e e e 


: 5 . 38 ) 


Substituting (5.20) in (5.37) one obtains 


®Cj ^Ci,Cj ^Ci + ’^Ci,Nk ®0i 


(5.39) 
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Similarly using (5.22) in (5.38), is given ly 


3 * f* 1 
e ee e 


m 

+ I 

k=l 


^e,Hk 


(5.40) 


Ihe second terms in (5.39) and (5.40) represent the 
net power flow into the system through the nodes of 
roup N and the equivalent respectively while the first 
terms give the net powers used up at the nodes of C and 
the equivalent respectively. Since the transfer admittance 
between nodes ’e’ and ’Nk’ in the equivalent is known from 
(5.35), by equating P^ and P^ and simplifying Y is 
obtained as 



T y s *. t* ^ . 

iii ik Cl .Cl, Co 




k=l i=l ^Ci,Nk 


K - - 

^ ®Ciki,lilk®nk^ • 


(5.41) 


Thus. (5.35) and (5.41) yield valu.es of the unknown 
parameters in the admittance matrix of (5*18). The elements 
in of (5.17) and (5.18) are unaffected. Prom these 
equations it is clear that for any coherent group, one 
row of elements has to be evaluated corresponding to 
node ' e* in the equivalent and involves only multi- 
plications, while the dynamic equivalent developed in 
reference [53] using the node elimination formula involves 
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matrix inversions. Considerable reduction in time is thus 
effected in the present method of equivalencing, 

5 . 4 NUMERIC A1 EXAMPLE: 

■The test system chosen for study is the 44 bus, 

15 machine and 56 lines system .sho-wn in Pigure 5.4. The 
data for this system is given an Appendix S. The region 
encompassing the buses 1,2,3»4 and 6 was assumed to be 
the study area. Buses 1 and 2 are generator buses while the 
rest are non-generator buses. The performance of the machines' 
1 and 2 in the study area was ana.lyzed when a fault on 
bus 6 was cleared after 0.1 second. 

In order to identify coherent machines, a detailed 
stability study was carried out on the unreduced system 
xfith the generator represented by a constant voltage behind 
their respective transient reactances. Since the purpose 
of this study is (i) to demonstrate the effectiveness of 
the coherency criterion developed in Section 5.2 and 
(ii) the feasibility of the coherency based dynamic equiva- 
lent derived in Section 5.3, to large scale power systems 
the computations of the quantities on the left hand sides 
of equations (5.14) for coherency recognition at the end 
of every time step was included in the program. The output 
lists out pairs of generators that are coherent according 
to the tolerences specified. Various values of and ^2 




FIG. 5.4: 44 BUS SAMPLE SYSTEM 
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wore tried and their effect on cohorencj?' was investigated. 
Table 5.1 gives the coherent groups for different tolerances 
assigned to and £ 2 * With a value of 0.0001 assigned to 
each, only one group consisting of machines 15,14,15 in 
the external system was found as coherent. This implies 
that these machines are tightly coupled. Swing curves 
plotted for these machines in Figure 5.5 show perfect 
coherency between them. A relaxation in £^ or £2 yields 
another coherent group comprising machines 9 and 12. 

Swing curves for these machines are shown in Figure 5.6. 

A rough visual inspection of the swing curves of these 
machine groups (13,14,15) and (9,12) would indicate that 
all these machines are coherent. However, this is 
critically dependent upon the values of and 
(cases 2 and 5 in Table 5.1). With £^ and £2 as above 
two distinct groups are indica,ted. Loosening the res- 
triction by allowing £-, and £., 0.0005 makes all 5 machines 

coherent. This study thus reveals the importance of a 
proper choice of ^1 .and £2 for detecting coherency. 

Having obtained the coherent groups, the data 
pertaining to these machines is fed to the equivalencing 
program that yields the djmamic equivalent of the system. 
The program computes the voltages, inertias, power input 
and self and transfer admittances of the equivalent 
generator replacing the coherent group. For this, the 
two groups obtained in case 2 of Table 5.1 were used. 
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■Tatle 5.1: Effect of and ^2 
coherency. 


Case No, 

"l 

^2 

MCHINBS D 
COHERENT GROUl 
1 

_ 

o 

1 



0.0001 

0.0001 

13,14,15 

— 

2 

0.0001 

0.0005 

^ 13,14,15 

9,12 

5 1 

1 

0.0005 

1 

0.0001 1 

i 

13,14,15 

9,12 

1 

4 

i 1 

0.0005 

i 

0.0005 1 
i 

9,12,13,14,15 
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A stability study was again conducted with these coherent 
groups replaced by their eQ’a.iva,lent in order to obtain the 
behaviour of machines 1 and 2 in the study area. Figures 
5.7 and 5.8 show the swing curves of these two machines 
obtained from the detailed study and using the equivalents. 
Close agreement is observed. The swing curves of the 
equivalent machines replacing the coherent gi-oups are 
shown in Figures 5'. 5 and 5.6, These curves indicate that 
the equivalent reacts faster than the machines it replaces. This 
is due to the fact that the electrical distance of this 
equivalent generator from the disturbance is reduced. But 
the nature of the swing is similar to the machines in the 
corresponding coherent group. Since the behaviour of the 
machines 1 and 2 in the study area is of concern, the 
equivalent representing the machines in the external area 
can be considered to have behaved well. The swing curves 
of machine '5' retained in the external area is also 
shown in Figure 5.9 obtained from a’ stability study on 
the reduced and unreduced .system. G-ood agreement is 
observed. 

5.5 COUGLUSIOUS; 

The necessity of an accurate equivalent that 
simulates the behaviour of the original system for 
stability studies has been felt in the recent years. An 
effective method to detect coherency in the system based 
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on energetical considerations has been presented in 
this chapter. The responses of the machines in a coherent 
group obtained by this technique, wore found to be 
perfectly matching as is evident from the swing curves. 

The dynamic equivalent derived utilizing the property 
of coherency is found to yield responses for machines in 
the study area, identical to those from a detailed 
simulation. Such an equivalent is extremely useful in 
the stability analysis of large scale power systems 
using Lyapunov methods since the number of unstable 
equilibrium points and hence their computation is 
greatly reduced. 



OHAPTHIR 


COxTOLUSIOITS 

6,1 STOIMRY AM COSTCIUSlOIjrxS 

The development of Lyapunov theory for the stability 
analysis ox power systems has been carried out with sufficient 
interest in recent years. The method has proven to be an 
attractive alternative to the simulation techniques both from 
the theoretical and practical points of view. T//hile sufficient 
work has been carried out for the construction of Lyapunov 
functions, the application of the method to realistic power 
systems from the practical point of view has been slow. Apart 
from the general conservative nature of the method inherent 
in the Lyapunov theorem itself, there are still certain 
drawbacks still to be overcome. One of these drawbacks is 
in regard to multimachine system modelling. It has been 
difficult to include effects like saturation flux decay, 
governor controls and transfer conductances in constructing 
a Lyapunov function. Hence the models were restricted to 
the simplest type. Therefore, alternate approaches are 
called for in order to consider these factors in the 
construction of a suitable Lyapunov function. Some of 
these include the use of vector Lyapunov fxinctions and 
the use of 'simplified models called the dynamic equivalents. 



Anotlier drawback of the Lyapunov method for the 
stability analysis of power systems has been in regard 
to the determination of the sta.bility boundary. The 
problem involves the computation of a large number of 
unstable equilibrium points and demands a considerable 
amount of computer time. This difficulty has been 
obviated to a large extent in a recent work of Prabhakara 
and El-Abiad. Improvements on this method in obtaining 
better stability boundaries are possible. 

In this research, emphasis has been, therefore, given 
to some aspects of the above two problems. The following 
are some of the conclusions arrived at. 

The nature of power system state models for the 
construction of lyaptmov functions has been a topic of 
considerable controversy over the past few years. It 
has been agreed that some of these models are uncontrol- 
lable and unobservable and hence cannot be used for the 
construction of Lyapunov functions using systematic pro- 
cedures of Kalman or the Moor e-Anderson theorem. Although 
this controversy is by now settled via the control 
theoretic concepts of minimal realization, reasons for 
uncontrollability in some of the models have not fully 
been explained. It has been shown in Chapter II that 
there is a physical dependency existing among the state 
variables in such models resulting in overdescription 
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and this indeed seems to be the cause of the so-called 
un.controllability of some of the models used in the 
literature. 

The new approach of "vector lyapunov function 
introduced in Chapter III is a two level concept and 
utilizes the decomposition and aggregation techniques 
for the stability analysis of large scale power systems. 
.Decomposition of the system is a key step in the appli- 
cation of these concepts. A procedure for decomposing 
the power system into lovr order subsystems has been 
proposed. The nature of decomposition is such that each 
of these oohsists of two machine groups i. It is therefoire 
possible to include the system refinements in each of 
these subsystems at the lower level and is a topic of 
further investigation. Tor the simple model considered 
at the subsystem level, lyapunov f"unctions hp.ve been 
constructed. Using the properties of these lyapunov 
functions and the na.ture of interactions among the sub- 
systems, the application of vector Lyapunov function at 
the higher hierarchiel level has been demonstrated. 
Preliminary results indicate that the method yields more 
conservative stability regions. However, the possibility 
of including system refinements outweighs this disadvantage. 
Also recent advances in large scale system stability theory 
indicate that sharper results may be obtained. 
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The subject roatter of Chapter IV has dealt with 
the determination of the stability boundary of a power 
system by a- method that s-voids the computation of the 
unstable eq_uilibrium points. The use of the sector 
violation information of the nonlinearities in the system 
model (which is in the Lure'-— Popov form) in obtaining the 
stability region for a multinonlinear system has been 
discussed. The results obtained are shown to be a 
generalization of the work of a similar nature in 
connection with the estimation of stability domain for a 
system with a single nonlinearity. On the basis of these 
results, the problem of estimating the stability region for 
the jaulti-machine power system is analyzed, Numerical 
computations on the 5-machine example presented and on 
other systems [79] indicate that the method, yields 
conservative estimates of the domains. But the simpli- 
city of computation is more attractive. 

Chapter V has presented a new method of identi- 
fying coherency and developing dynamic equivalents 
using energy functions. It frequently happens that 
certain section say ’A’ of the system is not adversely 
! affected by a disturbance in another section *B'. In 
such a case some of the machines in' A' tend to swing 
almost together during the period under study in groups. 
Such groups of coherent machines can be replaced by a 
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single equivalent machines for each group. A reduced 

f f 

ox'der dynamic equivalent of the system A can then he 
derived. In this chapter^ coherency identification and 
the development of the equivalent are based on certain 
energy criteria. Por this purpose an energy type lypunov 
function is used. An advantage of the method is that the 
angular differences at any time during the period of study 
is always compared with those of the initial angular 
differences. Results of tests conducted on the 44 bus 
system show that the method is extremely effective in 
respect of coherency recognition. The dynamic equivalent 
produces swing curves for the retained machines in almost 
full agreement with those obtained by a stability study 
of the original system, 

6.2 SU&GBSTIORS FOR FUTURB P^SEARCH: 

On the basis of the investigations carried out 
in this thesis, some of the unsolved problems that need 
further attention are briefly reviewed now. These are: 

1. The decomposition achieved in Chapter III is 
purely mathematical. Bach of the subsystems so obtained 
represents the dynamics of a 2-machine group. The number 
of subsystems derived thus becomes large as the number of 
machines increase. Alternative methods of decomposition, 
may be tried. 
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2, The stability region for the overall system is 
derived via a scalar Lyapimov function.. An alternative 
method of obtaining this region employing the stability 
boimdaries of the individual subsystems, as suggested by 
Weissenberger [ 77 ] may be persued. 

3, The conventional model of the system with a constant 
voltage behind the transient reactance has been only used 
in Chapter III. This simplicity in, the model has been 
retained in order to dojonstrate the applicability of the 
method to the power system problem. Further work is needed 
to include refinements like the governor control, voltage 
regulators, saliency etc., The method of decomposition 
proposed in this chapter will be useful for this purpose, 

4, The concept of a corner point in Chapter IT yields 
better stability domains than those obtained by other 
methods in the litera.ture. It has been shown that this 
point can be identified with the machine going out of 
step with respect to all other machines. This, however, 
requires further theoretical investigation, A recent 
paper of Gupta and El-ilbiad [82] discusses improvement in 
the stability region by considering 2-macliines going out 
of step etc. Examination of this heuristic concepts in 
the light of the theoretical results of Chapter IT is a 
subject of further research. Also the concepts of open 
Lyapunov surfaces introduced by Willems may be tried to 
improve the stability boundary . 
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APPBHDn: A 

SOME DBEINITIOFS Ilf COlBTECTIOlf flllH THBOEEM 3.1 

In this appendix, some definitions of terms used 
in theorem 3.1 nre given [73]. 

DEFIITITIOIT A.li A real valued function* (r) belongs to 
class K if it is defined, continuous and strictly increasing 
on 0 <_ r <. H vanishes at r = 0. 

DSP Ilf IT ION A. 2: A function V(X) is called positive definite 
in a domain if there exists a real valued function ij)(r) 
"belonging to class k such that 

( II 2t1 ) ^ 7(X) , X (A.l) 

DBEIhlTION A. 3: A fimction V(l) is called decrescent if 
there exists a function ^(r) of class k such that in the 
neighbourhood of the origin and for all t A t^ 

7(1) < ( 11 X II ) (A. 2) 

DEPIMITION A.4j A function 7(X) is said to be negative 
definite if -7(X) is positive definite. 

DEFIIITIOM A, 5- An autonomous system 

X = E(X) ; E(0) = 0 (A. 3) 

is said to be 'stable’ in the sense of Lyapunov if for 

each e >0, there exists an r| > 0 such that 
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II X(t) [| < s (A. 4) 

is valid whenever 

l!x(t^)|< n, .. 

where X(t) is a solution of (A. 5). 

DEfIKITION A, 6% The autonomous system (A. 3) is 
asymptotically stable whenever the origin X = 0 is 
stable and furthermore for some e > 0, every solution 

X 0 S-S t-><» . 
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ilPPEMDIX B 

PROOF OP TH30BEM 3.2 OP CHAPTER III 


In this appendix, the proof of theorem 3.2 is 
given [48]. 

Let us introduce a decrescent positive definite 
function V (V ) s 

S V 

(B.l) 

where h is a constant s-vector (h^, . • . h^) and is 

positive element by element which we denote by b > 0. 
Consider ^ S'S a candidate for the Lyapunov’s 
fimction for the composite system 5 . Taking the total 
time derivative along the solution X(t ; t^, of S we 

get 

vyip = h ■ (B.2) 

T 

Premultiplying the inequality (3.14) by b > 0 we get 

Vg(l^) < b^ R I . (B.3) 

Since R has non-negative off diagonal elements and 
satisfies conditions (3.18) R is Hurwitz [74]. Therefore 
conditions (3.18) are necessary and sufficient for the 
existence of an s-vector b > 0 for any s-vector c_> 0 
9uch that 


(B.4) 
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This implies 

^g(V < I (B.5) 

Since c. > 0 and ¥ >0, <0. Also wSince _b > 0, 

^ Thus the glohal asymptotic sta,bility 

S •“'V' 

of S is assured. This completes the proof of the 
theorem 3.2 of Chapter III. 
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APPENDIX C 

res'ljlis op appeigation op the pioorb-;\itdehson' s theorem 

TO THE PO¥ER SYSPBM PROBLEM 


In this appendix the theoren due to Moore-Anderson 
[8l] is stated without proof. Results of application of 
this theorem to the power system in a manner outlined 
in reference [10 ], to derive a lya,punov function '7(Z) 
of the quadratic in state variables plus integral of the 
nonlinearity type are included. 


THE MOORE-AKDERSOIT THEOREM: If there exist diagonal matrices 
a = diag(aj^), i = 1,2,... m and Q = diagCg^^), i = l,2,,..m 

with a - >. 0, q . >. 0 and a . +n . ^ 0 and -ix./a. is not a 
1 1-^1 1 

pole of the ith row of ¥(s) such that 

Z(s) = a 0“^ + (a + Q(s)) ¥(s) (C.l) 


is positive real, then the system described by the set 
of equations 


X = A X ~ B f (o) 
n = X 


(C.2) 


is stable. 

In equation (4.29), the matrix G- = diag(gj^) , 
1 — “ 1 , 2 , . . . 


m with 
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0 < f^(cf^) < gj_ cy2 for gj_> 0, i = l,2,...in. 

(C.3) 

For the application of the theorem one needs the 
following definition of positiTe real f -unctions; 
DEFINITION Q.l; An mxn matrix Z(s) is positive real if 
and only if the following conditions are satisfied. 

1. Z(s) is rea.1 rational 

2. Z(s) has no poles in Re s> 0 

3. Poles of Z(s) on Re s = 0 are simple 

4. For each pole on Re s = 0, the residue matrix R 
is hermit ian and positive semidef init e, 

and 

5. The hernitian pa,rt of Z(s) is positive 

somidef inite, i.e. 

Zg(jco) + Zg(~3aj) > 0 for all real (a 


Proof of the above theorem is contained in reference [81], 


Utilizing theorem 1 a.nd the above definition on 
positive real functions systematic construction of 
lyapunov function for the system (4.40) resulted [lO] 
in a lyapunov function of the form (4.41) ■'/hiere the matrix 
P given by 




(C.4) 


and are 


i 
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(i) Uniform damping cases 

'‘ip * 

7® ! 
(n-l)(ii-.i)| 

where 


(0.5) 


^2 I M . -m.M^ 

j=l ^ - 2 3 


'(n-l)(n-l) 


.1 M-f 


“M^M2 


"3,1 '■>3 


•• -“2«n 


. . -M^M 
3 n 


-M M 
n 2 


n 3 


•• ""n I 
a=i 


(n-l) (n-1) 


a f 

^ (n-l) (n-l) 


( 0 . 6 ) 

(0,7) 


(n-l) (n-l) ^'^( n-l) (n-l) 


( 0 . 8 ) 


T- for all i = l,2,...m 


(ii) Fon-nniform damping 


case; 


Here 


(0.9) 


(2n-l)(2n-l) 


n(n-l) 


n(n-l) I Xn-l)(n-l) 


( 0 . 10 ) 
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where 


n 






I D 

i=l ^ 


^2 I ®i 
0=1 ^ 
jj^2 


"^2^3 


•D 2 ^Df^ 


n 

j=l 

oVl 


-V2 


■D2^n 




-D I™ 
n 3 


n^l 

(C.ll) 


^ n(n--l) ^ rm ^^n(rL-l) (n-l) (n-1) 


n 


and ^ = ( I -^) M . 

nn - ^-1 A. 


nn 


. X. 

1=1 1 


( 0 . 12 ) 

(0.13) 


The constant is given hy 


n 


Pi = 2( I T- ) 

^ i=l i 


(0.14) 
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APPJJJSTDIZ D 

DATA POH THE 5 ffilCHIHS .SYSTEM 

The line constants and ras-chine do,ta are given in 
Tables D.l and D.2. The load flow data prior to the 
fault on the • .-system is given in Table S.3. All data 
are given in p.u. 


Table D.l: Line Data 


Line 

Ho. 

Bus code 

from bus to bus 

Imped 

H 

anc e 

X 

Shunt 

susceptanc 

B 

1 

1 

2 

0.0092 

0.0582 

0.0348 

2 

2 

3 

0.0000 

0.0400 

0.0000 

3 

3 

4 

0.0142 

0.0970 

0.3490 

4 

2 

8 

0.0038 

0.0408 

0.0090 

5 

2 

7 

0.0262 

0.1094 

0.0060 

6 

2 

6 

0.0204 

0.1420 

0.0270 

7 

2 

5 

0.0310 

0.0882 

0.0050 

8 

5 

6 

0.0200 

0.1236 

0.0080 

9 

6 

7 

0.0118 

0.0643 

0.0227 

10‘ 

8 

9 

0.0080 

0.0180 

0.0900 

11 

8 

10 

0.0350 

0.1340 

0.0260 

12 

10 

11 

0.0350 

0.1340 

0.0265 

13 

9 

11 

0.0010 

0.0040 

0.0170 

14 

10 

12 

0.0196 

0.0610 

0.1260 

15 

1 

12 

0.0240 

0.0862 

0.0203 
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Table D.Zi 

Ma.chine late. 


Machine 

lo. 

At bus 

Inertia 

constant 

(H) 

Squix’-alent 
( p.u. 
R 

impedance 

) 

Z 

1 

1 

2.30 

0.0420 

0.1570 

2 

4 

o 

• 

O 

H 

0.0000 

0.1361 

3 

6 

12.72 

0.0000 

0.0717 

4 

9 

27.00 

0.0000 

0.0326 

5 

12 

00 

0.0060 

0.0280 


Table D.3* Pre-fault load Flow Conditions. 


Bus voltage G-enero-tion Load 


Fo,. 

Magnitude 

Angle 
(deg. ) 

Heal 

Reactive 

Real 

Reactive 

1 

0.996 

0.1 

0.900 

0.091 

0.170 

0.000 

2 

1.007 

-2.0 

0.000 

0.000 

1.400 

0.200 

3 

0.912 

2.8 

0.000 

0.000 

0.000 

0.000 

4 

1.000 

12.8 

1.760 

0.442 

0.000 

0.000 

5 

1.003 

-3^8 

0.000 

0.000 

0.900 

0.000 

6 

1.024 

0.1 

3.200 

0.511 

1.260 

0.150 

7 

0.999 

-3.7 

0.000 

0.000 

1.350 

0.200 

8 

0.992 

-3.2 

0.000 

0.000 

0.990 

0.400 

9 

1.000 

-2.9 

3.500 

1.462 

2.500 

0.900 

10 

0.963 

-2.5 

0.000 

0,000 

0.560 

0.200 

11 

0.998 

-2..0 

0.000 

0.000 

0.560 

0.200 

12 

0.956 

0.0 

0,443 

-0.723 

0.000 

0.000 
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APPENDIX E 

DATA POR THE 44 BUS SYSTEM 

The bus conditions in the pre-fault state are given 
in Table B.l. Line and shunt capacitor data are shown 
in Tabl® B.2 and E,3 respectively. The generator data 
is given in Table E.4. All data is in p.u, on a 200 MVA 
base. Damping constants of all machines are neglected. 

Table E.l: Pre-fault load Plow 
Conditions . 


Bus Bus voltage Generation load 


No. 

Magnitu- 

de 

Phas e 
angle 
(deg. ) 

Real 

Reactive Real 

Reactive 

1 

1.0300 

0.0000 

4.1087 

0.4539 

0.0000 

0.0000 

i 

1.0757 

-4. 2430 

0.0000 

0.0000 

0.0000 

0.0000 

3 

1.0250 

-2.4140 

2.3000 

0.8542 

0.0000 

0.0000 

4 

1.0503 

-6.6600 

0.0000 

0.0000 

0.0000 

0.0000 

5 

1.0497 

-12.8020 

0.0000 

0 . 0000 

0.8500 

0 . 6000 

6 

1.0245 

il6.8260 

0.0000 

0.0000 

0 . 6000 

0.3700 

7 

1.0000 

-3.0020 

1.5000 

0.5200 

0.0000 

0 .0000 

8 

1.0228 

-8.0380 

0.0000 

0.0000 

1.5400 

0.7450 

9 

1.0250 

-9.6280 

0.0000 

0.0000 

0.2683 

0 . 2010 

10 

1.0378 

-14.3620 

0.0000 

0.0000 

1.5000 

1.0005 

11 

1.1892 

-18.3320 

0.0690 

0.0000 

0.0000 

0.0000 

12 

1.1348 

-21.9430 

0.0000 

0.0000 

1.1750 

0.8000 

13 

1.0240 

- 1.4560 

0.4950 

0.1793 

0.0000 

0.0000 

14 

1.0405 

-7.6740 

0.0000 

0 . 0000 

0.0900 

0.0435 

15 

1.0390 

-33.3520 

0.1000 

0.0620 

0.0000 

0.0000 

16 

1.0800 

-19.4570 

1.0120 

0.4096 

0.0000 0.0000 
( continued) 
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Taole E.l continued. 


17 

1.0529 

- 23.4400 

0.0000 

0.0000 

1.2500 

0.7500 

18 

1.0594 

- 22.54 60 

0.0000 

0.0000 

0 . 9000 

0.6500 

19 

1.0297 

- 27.3030 

0.0000 

0.0000 

0 . 2000 

0.1500 

20 

1.0354 

- 27.1950 

0.0000 

0.0000 

0.6255 

0.4690 

21 

1.0223 

- 36.4070 

0.0000 

0.0000 

0.2250 

0.0750 

22 

1.0077 

- 26.9160 

0.0000 

0.0000 

0.4500 

0.3380 

23 

1.0832 

-a .8800 

0.1225 

0.0000 

0.0000 

0.0000 

24 

1.0500 

- 12.5550 

0.8350 

0.1471 

0.0000 

0.0000 

25 

1.0646 

- 17.5740 

0.0000 

0.0000 

0.0000 

0 . oooo 

26 

1.0360 

- 20.3840 

0.0000 

0.0000 

0.2315 

0.1735 

27 

1.0098 

- 25.0910 

0.0000 

0.0000 

0.7570 

0.5650 

28 

0..9900 

- 20.1920 

1.2400 

0.8687 

0.0000 

0.0000 

29 

1.0154 

- 26.6640 

0 .0000 

0.0000 

0.9460 

0.7090 

30 

0.9890 

- 32.3770 

0.0000 

0.0000 

0.0000 

0.0000 

31 

1.0630 

- 7.2760 

0 . 5080 

0.1149 

0.0000 

0.0000 

32 

0.9926 

- 20.4150 

0.0000 

0.0000 

0.2750 

0 . 2100 

33 

1.1043 

- 12.1480 

0.0850 

0.0000 

0.0000 

0.0000 

34 

1.0462 

- 19 . 0310 

0.0000 

0.0000 

0.1000 

0.0750 

35 

1.0739 

- 12.9000 

0.0000 

0.0000 

0.2500 

0.1500 

36 

1.0723 

'- 16 . 5140 

0,0000 

0.0000 

0.4500 

0.2200 

37 

1.0450 

- 7.9030 

1.5000 

0.5398 

0.0000 

0.0000 

38 

1.0 673 

- 13.0120 

0.0000 

0.0000 

0.0000 

0.0000 

39 

0.9746 

- 24.4070 

0.0000 

0.0000 

. 0.3375 

0.2525 

40 

• 1.0455 

- 15. 1690 

0.0000 

0.0000 

0.0000 

0.0000 

41 

1.0390 

- 10.3350 

0.7350 

0.3650 

0.0000 

0.0000 

42 

1.0680 

- 13.2620 

0.0000 

0.0000 

0.0000 

0.0000 

43 

1.0400 

- 6.9430 

1.0800 

0.1122 

0.0000 

0.0000 

44 

0 . 993 ffi 

- 28.8730 

0.0000 

0.0000 

2.1975 

1.3490 
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Table 

B.2; line 

Data 



SI. 

Mo, 

3us God e 

fr.om btis to bus 

Impedan ce 

Res is- Reactance 

Off nominal 

1 -Q turns 
ratio 

1 

8 

7 

0.0000 

0.0570 

0.0000 

1.050 

2 

8 

14 

0.0220 

0.0520 

0.0135 

1.000 

3 

8 

9 

0.0520 

0.1270 

0.0140 

1.000 

4 

9 

10 

0.0660 

0.1610 

0.0180 

1.000 

5 

14 

9 

0.0270 

0.0700 

0.0070 

1.000 

6 

10 

12 

0.1570 

0.3860 

0.0630 

1.000 

7 

11 

12 

0.0000 

1.1730 

0.0000 

1.050 

8 

14 

13 

0.0000 

0.2220 

0.0000 

1.050 

9 

20 

17 

0.0250 

0.1250 

0.1130 

1.000 

10 

7 4 

3 

0.0000 

0.0330 

0.0000 

1.050 

11 

14 

4 

0.0000 

0.0530 

0.0000 

1.000 ■ 

12 

4 

10 

0.1600 

0.1310 

0.0710 

0.900 

13 

2 

4 

0.0000 

0.0620 

0.0000 

1.000 

14 

4 

6 

0.0580 

0.2860 

0.0645 

1.000 

15 

2 

1 

0.0000 

0.0190 

0.0000 

1.050 

16 

6 

17 

0.0340 

0.1670 

0.1500 

1.000 

17 

5 

12 

0.0440 

0.2670 

0.1010 

0.850 

18 

2 

5 

0.0050 

0.0510 

0.6706 

1.000 

19 

17 

16 

0.0000 

0.0780 

0.0000 

1.000 

20 

17 

18 

0.0130 

0.0640 

0.0580 

1.000 

21 

18 

5 

0.0620 

0.1232 

0.2012 

1.150 

22 

19 

22 

0.0840 

0.1880 

0.0210 

1.000 

23 

23 

22 

0.0770 

0.7610 

0.0215 

1.050 

24 

22 

27 

0.1090 

0.1960 

0.0220 

1.000 

25 

30 

27 

0.0000 

0.0800 

0.0000 

0.950 

26 

20 

29 

0.0210 

0.1030 

0.0920 

1.000 

27 

29 

28 

0.0000 

0.0830 

0.0000 

1.100 

28 

29 

44 

0.1700 

0,8400 

0.0760 

1.000 


( continusd) 
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29 

30 

29 

0.0370 

0.1950 

0.0390 

1.000 

30 

44 

39 

0.0290 

0.1520 

0.0300 

1.000 

31 

39 

32 

0.0160 

0.0850 

0.0170 

1,000 

32 

32 

34 

0.0560 

0.3040 

0.0035 

0.950 

33 

24 

25 

0.0000 

0.1200 

0.0000 

0.975 

34 

25 

26 

0.0370 

0.0900 

0.0100 

1.000 

35 

26 

27 

0.0830 

0.1540 

0.0170 

1.000 

36 

26 

34 

0.1070 

0.1970 

0 , 2100 

1.000 

37 

40 

36 

0.0000 

0.0800 

0.0000 

0.950 

38 

34 

36 

0.0900 

0 . 2310 

0.0060 

1.000 

39 

40 

43 

0.0230 

0.1420 

0.1040 

1.050 

40 

35 

31 

0.0000 

0.1880 

0.0000 

1.025 

41 

38 

37 

0.0000 

0.0630 

0.0000 

1.050 

42 

38 

32 

0.0280 

0.1440 

0.0290 

1,000 

43 

38 

42 

0.0080 

0.0420 

0.0080 

1.000 

44 

38 

40 

0.0250 

0.1190 

0.0240 

1.000 

45 

38 

35 

0.0000 

0.1600 

0.0000 

1.000 

46 

36 

35 

0.0740 

0.1880 

0.0190 

1.000 

47 

40 

44 

0.0074 

0.1974 

0.2415 

1.000 

48 

25 

22 

0.3490 

0.8900 

0.0120 

1.000 

49 

18 

19 

0.0800 

0.4650 

0.0420 

0,975 

50 

15 

21 

0 . 2900 

0.7830 

0.0025 

0.9500 

51 

17 

21 

0.3980 

1.8960 

0.0265 

0.975 

52 

41 

42 

0.0000 

0.0810 

0.0000 

0.950 

53 

42 

40 

0.0100 

0.0480 

0.0430 

1.000 

54 

40 

30 

0.0230 

0.2300 

0.0695 

1.000 

55 

34 

33 

0.1170 

1.7150 

0 .0010 

0.950 

56 

4 

6 

0.0580 

0.2860 

0,0645 

1.000 
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'Table B.3» Shunt Capacitor Data 


SI. Bo. 

Bus Bo. 

Admittance 

G 

B 

1 

5 

0.0 

0.3400 

2 

12 

o 

« 

o 

0.6100 . 

3 

17 

0.0 

0.3750 

4 

18 

0.0 

0.2900 

5 

19 

O 

• 

o 

0.1100 

6 

20 

0.0 

0.4200 

7 

22 

0.0 

0.2200 

8 

27- 

O 

• 

O 

0.1900 

9 

44 

o 

• 

o 

1.0000 



Table 

B,42 Machine Data 


SI, Bo. 

At bus 

Inertia constant. 

Transient 

reactance 

1 

1 

43.7500 

0.0309 

2 

3 

15.4050 

0.0689 

3 

7 

6.6480 

0 . 2010 

/. 

11 

0.3450 

1.1733 

5 

13 

1.9230 

0.5278 

6 

15 

0.5700 

2.0553 

7 

16 

8.2500 

0.1965 

8 

23 

0 . 6900 

1.6232 

9 

24 

2.2600 

0.3520 

10 

28 

4.2790 

0.1185 

11 

31 

2.9520 

0.4045 

12 

33 

0.4350 

2.0000 

13 

37 

6.1900 

0.1500 

14 

41 

4.8640 

0.2555 

15 

43 

H 

O 

0.2319 
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